Blunt-nosed Leopard Lizard Data Analysis: Behavior and Thermoregulation vs Water Physiology

Savannah Weaver

Packages

`%nin%` = Negate(`%in%`)
if (!require("tidyverse")) install.packages("tidyverse") 
library("tidyverse")# for dplyr workflow of calculations
if (!require("lmerTest")) install.packages("lmerTest") 
library("lmerTest") # LMMs
if (!require("ggplot2")) install.packages("ggplot2")
library("ggplot2") # plots
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # pretty colors
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("Hmisc")) install.packages("Hmisc")
library("Hmisc") # correlation matrix
if (!require("emmeans")) install.packages("emmeans")
library("emmeans") # marginal means & confints
if (!require("car")) install.packages("car")
library("car") # Anova
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown format

Set Colors & Themes

savs_ggplot_theme <- theme_classic() + 
                     theme(text = element_text(color = "black", 
                                               family = "sans", 
                                               size = 12),
                           axis.text = element_text(color = "black", 
                                                    family = "sans", 
                                                    size = 8),
                           legend.text = element_text(color = "black", 
                                                      family = "sans", 
                                                      size = 8),
                           legend.text.align = 0,
                           legend.position = "bottom"
                           )

savs_ggplot_theme2 <- theme_classic() + 
                     theme(text = element_text(color = "black", 
                                               family = "sans", 
                                               size = 10),
                           axis.text = element_text(color = "black", 
                                                    family = "sans", 
                                                    size = 10),
                           legend.text = element_text(color = "black", 
                                                      family = "sans", 
                                                      size = 8),
                           legend.text.align = 0,
                           legend.position = "bottom"
                           )

my_brew_cols_4 <- brewer.pal(11, "RdYlBu")[c(5,4,10,11)] # this one!

Background

This data was collected on Blunt-nosed Leopard Lizards (Gambelia sila) during their active season in 2021. In this document, we will compute metrics of thermoregulation and behavior, then assess whether differences in those variables could be related to differences in hydric physiology.

All of the data wrangling and preparation for the thermal and microhabitat use data is done within this document, but is done in the ‘wrangling_’ documents for the hydric data.

First, we compute maximum temperatures and thermoregulatory accuracy (thermal metrics) and microhabitat use for each lizard at different time intervals, then we pair that with our hydric physiology data, and look for any relationships.

Load & Wrangle Data

Lizard Data

This is the same dataframe used in ‘analysis_1’. For info on data preparation, see the ‘wrangling_’ files, and for an explanation of each of the variable names, see that section in ‘analysis_1’.

# measurement dates
apr_may_dates <-  as.Date(c("2021-04-23", "2021-04-24", "2021-04-25"))
may_dates <- as.Date(c("2021-05-07", "2021-05-08", "2021-05-09"))
  
liz_dat <- read_rds("./data/G_sila_clean_full_dat.RDS") %>%
  # we only use data for radio collared lizards in this analysis
  dplyr::filter(complete.cases(radio_collar_serial)) %>%
  # add date chunks to help look at values over time
  mutate(date_chunks = factor(case_when(capture_date %in% apr_may_dates ~ "Apr 26 - May 6",
                                 capture_date %in% may_dates ~ "May 10 - 20")),
         radio_collar_serial = factor(radio_collar_serial))
summary(liz_dat)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 03:39:00.00   Min.   :2021-04-23   229-069: 3         
##  1st Qu.:2021-04-23 07:54:30.00   1st Qu.:2021-04-23   245-742: 3         
##  Median :2021-04-24 03:50:30.00   Median :2021-04-25   245-744: 3         
##  Mean   :2021-04-28 07:03:08.27   Mean   :2021-05-12   245-745: 3         
##  3rd Qu.:2021-05-07 04:41:30.00   3rd Qu.:2021-05-08   245-748: 3         
##  Max.   :2021-05-08 07:19:00.00   Max.   :2021-07-14   252-883: 3         
##  NA's   :18                                            (Other):58         
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:76          F-12   : 3    Female:32   Min.   :25.00   Min.   : 85.0  
##  Class :character   M-09   : 3    Male  :44   1st Qu.:33.70   1st Qu.: 96.0  
##  Mode  :character   M-10   : 3                Median :38.00   Median :100.0  
##                     M-11   : 3                Mean   :39.20   Mean   :101.7  
##                     M-19   : 3                3rd Qu.:43.45   3rd Qu.:108.0  
##                     M-20   : 3                Max.   :58.80   Max.   :122.0  
##                     (Other):58                                               
##         tmt     tmt_mass_change_g capture_to_msmt  hematocrit_percent
##  Water Tmt:38   Min.   :-3.0000   Min.   :  11.0   Min.   :17.00     
##  Sham Tmt :38   1st Qu.:-1.0000   1st Qu.:  62.5   1st Qu.:30.00     
##  No Tmt   : 0   Median : 0.0000   Median : 152.0   Median :34.00     
##                 Mean   :-0.2262   Mean   : 401.1   Mean   :34.35     
##                 3rd Qu.: 0.0000   3rd Qu.: 265.2   3rd Qu.:38.50     
##                 Max.   : 3.0000   Max.   :1665.0   Max.   :58.00     
##                 NA's   :15        NA's   :18       NA's   :1         
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :24.00   Min.   : 0.650   Min.   :19.07   Min.   :12.20  
##  1st Qu.:30.00   1st Qu.: 7.989   1st Qu.:28.02   1st Qu.:14.79  
##  Median :32.00   Median : 9.620   Median :29.31   Median :20.38  
##  Mean   :31.92   Mean   : 9.854   Mean   :29.18   Mean   :22.70  
##  3rd Qu.:34.00   3rd Qu.:12.324   3rd Qu.:31.74   3rd Qu.:27.49  
##  Max.   :36.00   Max.   :21.673   Max.   :33.55   Max.   :40.22  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.205   Min.   :1.415   Min.   :309.3           April:39   1 :39  
##  1st Qu.:3.784   1st Qu.:2.487   1st Qu.:345.0           May  :23   3 :23  
##  Median :4.076   Median :3.229   Median :361.0           July :14   12:14  
##  Mean   :4.100   Mean   :3.212   Mean   :363.6                             
##  3rd Qu.:4.685   3rd Qu.:3.982   3rd Qu.:377.9                             
##  Max.   :5.188   Max.   :4.319   Max.   :436.5                             
##                                  NA's   :1                                 
##     week_num      ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   : 1.000   Min.   :2021-04-24   Gravid    :16   Min.   :3.000  
##  1st Qu.: 1.000   1st Qu.:2021-04-24   Not Gravid:12   1st Qu.:3.750  
##  Median : 1.000   Median :2021-04-24   NA's      :48   Median :4.000  
##  Mean   : 3.632   Mean   :2021-04-29                   Mean   :3.812  
##  3rd Qu.: 3.000   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :12.000   Max.   :2021-05-08                   Max.   :5.000  
##                   NA's   :48                           NA's   :60     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage   
##  Min.   :0.5700          Min.   :1.710                small round: 1  
##  1st Qu.:0.8425          1st Qu.:2.385                large round: 9  
##  Median :0.9950          Median :2.885                soft       : 2  
##  Mean   :1.1131          Mean   :3.148                soft oblong: 1  
##  3rd Qu.:1.4200          3rd Qu.:3.985                firm oblong: 2  
##  Max.   :1.8000          Max.   :4.740                hard oblong: 1  
##  NA's   :60              NA's   :60                   NA's       :60  
##    dev_point         SMI                date_chunks
##  Min.   :1.00   Min.   :23.67   Apr 26 - May 6:39  
##  1st Qu.:2.00   1st Qu.:32.29   May 10 - 20   :23  
##  Median :2.00   Median :34.40   NA's          :14  
##  Mean   :2.75   Mean   :34.77                      
##  3rd Qu.:3.25   3rd Qu.:37.43                      
##  Max.   :5.00   Max.   :45.63                      
##  NA's   :60
# subsets by date
liz_dat_apr_may_only <- liz_dat %>%
  dplyr::filter(complete.cases(date_chunks))
liz_dat_july_only <- liz_dat %>%
  dplyr::filter(month == "July")

# ID info only
liz_IDs_only <- liz_dat %>%
  group_by(radio_collar_serial, individual_ID, sex_M_F) %>%
  summarise(n = n()) %>%
  dplyr::select(-n)
## `summarise()` has grouped output by 'radio_collar_serial', 'individual_ID'. You
## can override using the `.groups` argument.

Validation

This is a dataframe I can use to make sure we do not have data for these lizards after these dates, because they were either found dead or their dropped collars were found.

check_gone <- data.frame(Serial = factor(c("229-072", # found dead June 7
                                           "229-060", # collar found May 18
                                           "229-074", # collar found April 28
                                           "229-070", # became -a on May 8
                                           "245-762", # became -a on May 8
                                           "229-073", # dropped sometime before May msmts
                                           "245-751"  # died May 7/8
                                           )),
                         date_def_gone = as.Date(c("2021-06-07", "2021-05-18",
                                                   "2021-04-28", "2021-05-07",
                                                   "2021-05-07", "2021-05-07",
                                                   "2021-05-07"),
                                                 format = "%Y-%m-%d"))
check_gone
##    Serial date_def_gone
## 1 229-072    2021-06-07
## 2 229-060    2021-05-18
## 3 229-074    2021-04-28
## 4 229-070    2021-05-07
## 5 245-762    2021-05-07
## 6 229-073    2021-05-07
## 7 245-751    2021-05-07

Tsurf Wrangling

We measured the surface body temperature of lizards throughout the study using Holohil temperature-sensitive transmitters, which is the data we load in this section.

variables in this dataframe:

  • DateTime = date and time the temperature was logged
  • Serial = the radio-collar serial number for that observation (to be used to link to specific individuals)
  • Freq = what radio frequency the collar was on
  • pp = the pulse period recorded by the radio-collar receiver (pulse period is a proxy for temperature)
  • ppTempC = the pp value converted into temperature in degrees Celsius using third order polynomial functions estimated from the reference transformations provided by Holohil
temps <- read.csv("./data/telemetry/tower_temps_all.csv") %>%
  # properly format data classes
  mutate(date_only = as.Date(substr(DateTime, 1, 10), format = "%Y-%m-%d"),
         hour_of_day = as.numeric(substr(DateTime, 12, 13)),
         DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M:%S"),
         Serial = (paste(substr(Serial, 1, 3), 
                               substr(Serial, 4, 7), sep = "-")),
         Freq = round(Freq, 2))

# look at temps dataframe
summary(temps)
##     DateTime                         Serial               Freq      
##  Min.   :2021-04-25 13:11:35.00   Length:107452      Min.   :150.0  
##  1st Qu.:2021-05-02 18:18:19.75   Class :character   1st Qu.:150.4  
##  Median :2021-05-11 07:58:47.00   Mode  :character   Median :151.1  
##  Mean   :2021-05-15 02:34:25.98                      Mean   :151.0  
##  3rd Qu.:2021-05-17 22:48:35.50                      3rd Qu.:151.5  
##  Max.   :2021-07-13 13:11:42.00                      Max.   :152.0  
##     ppTempC             pp         date_only           hour_of_day   
##  Min.   :-23.99   Min.   :1001   Min.   :2021-04-25   Min.   : 0.00  
##  1st Qu.: 31.80   1st Qu.:1565   1st Qu.:2021-05-02   1st Qu.:10.00  
##  Median : 36.26   Median :1703   Median :2021-05-11   Median :13.00  
##  Mean   : 35.20   Mean   :1734   Mean   :2021-05-14   Mean   :12.69  
##  3rd Qu.: 39.06   3rd Qu.:1885   3rd Qu.:2021-05-17   3rd Qu.:16.00  
##  Max.   : 76.93   Max.   :2799   Max.   :2021-07-13   Max.   :23.00
unique(temps$Serial)
##  [1] "229-059" "229-060" "229-063" "229-065" "229-066" "229-069" "229-070"
##  [8] "229-072" "229-073" "229-074" "229-077" "229-078" "229-079" "229-080"
## [15] "245-741" "245-742" "245-745" "245-748" "245-751" "245-752" "245-753"
## [22] "245-754" "245-756" "245-757" "245-759" "245-760" "245-761" "245-762"
## [29] "252-881" "252-882" "252-883" "252-884" "252-885" "252-886" "252-887"

Check whether we have observations on lizards that were not alive or were not actually wearing their radio-collars because they dropped them.

temp_check <- temps %>%
  left_join(check_gone, by = 'Serial') %>%
  dplyr::filter(date_only > date_def_gone) %>%
  arrange(Serial) %>%
  dplyr::select(Serial, date_only, date_def_gone)
unique(temp_check$Serial)
## [1] "229-060" "229-070" "229-072" "229-073" "229-074" "245-751" "245-762"

Great… we have data for all 7 individuals that we need to clean up for

temps2 <- temps %>%
  # delete the obs for all Serials but 229-070 & 245-762
  dplyr::filter(!(Serial == "229-072" & date_only >= "2021-06-07")) %>%
  dplyr::filter(!(Serial == "229-060" & date_only >= "2021-05-18")) %>%
  dplyr::filter(!(Serial == "229-074" & date_only >= "2021-04-28")) %>%
  dplyr::filter(!(Serial == "229-073" & date_only >= "2021-05-07")) %>%
  dplyr::filter(!(Serial == "245-751" & date_only >= "2021-05-07")) %>%
  # change the serial and ID of F-08 and M-03 to be "-a" bc collar was switched
          # for F-08a
  mutate(Serial = replace(Serial, # what I want to change
                          Serial == "229-070" & # if this
                          date_only > "2021-05-07", # and this
                          "229-070a"), # change to this
         # for M-03a
         Serial = replace(Serial, # what I want to change
                          Serial == "245-762" & # if this
                          date_only > "2021-05-07", # and this
                          "245-762a"), # change to this
         Serial = factor(Serial))

Check that it worked:

temps_doublecheck <- temps %>%
  # SAVE the obs for all Serials but 229-070 & 245-762
  dplyr::filter((Serial == "229-072" & date_only >= "2021-06-07") |
                  (Serial == "229-060" & date_only >= "2021-05-18") |
                  (Serial == "229-074" & date_only >= "2021-04-28") |
                  (Serial == "229-073" & date_only >= "2021-05-07") |
                  (Serial == "245-751" & date_only >= "2021-05-07"))
nrow(temps_doublecheck) + nrow(temps2) == nrow(temps)
## [1] TRUE
temp_check2 <- temps2 %>%
  left_join(check_gone, by = 'Serial') %>%
  dplyr::filter(date_only > date_def_gone) %>%
  arrange(Serial) %>%
  dplyr::select(Serial, date_only, date_def_gone)
temp_check2
## [1] Serial        date_only     date_def_gone
## <0 rows> (or 0-length row.names)
unique(temps2$Serial)
##  [1] 229-059  229-060  229-063  229-065  229-066  229-069  229-070a 229-070 
##  [9] 229-072  229-073  229-074  229-077  229-078  229-079  229-080  245-741 
## [17] 245-742  245-745  245-748  245-751  245-752  245-753  245-754  245-756 
## [25] 245-757  245-759  245-760  245-761  245-762a 245-762  252-881  252-882 
## [33] 252-883  252-884  252-885  252-886  252-887 
## 37 Levels: 229-059 229-060 229-063 229-065 229-066 229-069 229-070 ... 252-887

Subset the temperature data so that we do not use the temperatures from during hydric measurement periods. Also add columns for time interval chunks to make averaging easier later.

# capture weekend dates
cap_dates <- as.Date(c("2021-04-23", "2021-04-24", "2021-04-25", # April
                        "2021-05-07", "2021-05-08", "2021-05-09"), # May
                     format = "%Y-%m-%d")

# grouping notes by date
group_dates <- data.frame(date_only = c(seq(as.Date("2021-04-26"), 
                                            as.Date("2021-07-13"), by = 1))) %>%
  dplyr::filter(date_only %nin% cap_dates) %>%
  mutate(date_chunks = factor(c(rep("Apr 26 - May 6", 11), 
                                rep("May 10 - 20", 11), 
                                rep("May 21 - 31", 11),
                                rep("Jun 1 - 11", 11),
                                rep("Jun 12 - 22", 11),
                                rep("Jun 23 - Jul 3", 11), 
                                rep("Jul 4 - 13", 10)),
                              levels = c("Apr 26 - May 6", "May 10 - 20", 
                                         "May 21 - 31", "Jun 1 - 11", "Jun 12 - 22", 
                                         "Jun 23 - Jul 3", "Jul 4 - 13")))

# subset dataset and add grouping notes
temp_sub <- temps2 %>%
  dplyr::filter(date_only %nin% cap_dates) %>%
  left_join(group_dates, by = 'date_only')
summary(temp_sub)
##     DateTime                          Serial           Freq      
##  Min.   :2021-04-26 00:01:18.00   245-756: 5195   Min.   :150.0  
##  1st Qu.:2021-05-02 12:14:29.75   245-753: 5183   1st Qu.:150.4  
##  Median :2021-05-12 02:49:33.50   229-078: 4971   Median :151.1  
##  Mean   :2021-05-15 12:15:37.34   252-886: 4487   Mean   :151.0  
##  3rd Qu.:2021-05-18 14:59:10.50   229-069: 4457   3rd Qu.:151.5  
##  Max.   :2021-07-13 13:11:42.00   252-881: 4182   Max.   :152.0  
##                                   (Other):68779                  
##     ppTempC             pp         date_only           hour_of_day   
##  Min.   :-23.99   Min.   :1001   Min.   :2021-04-26   Min.   : 0.00  
##  1st Qu.: 32.02   1st Qu.:1561   1st Qu.:2021-05-02   1st Qu.:10.00  
##  Median : 36.40   Median :1698   Median :2021-05-12   Median :13.00  
##  Mean   : 35.33   Mean   :1729   Mean   :2021-05-14   Mean   :12.61  
##  3rd Qu.: 39.14   3rd Qu.:1877   3rd Qu.:2021-05-18   3rd Qu.:16.00  
##  Max.   : 76.93   Max.   :2798   Max.   :2021-07-13   Max.   :23.00  
##                                                                      
##          date_chunks   
##  Apr 26 - May 6:39496  
##  May 10 - 20   :33790  
##  May 21 - 31   : 9326  
##  Jun 1 - 11    : 6163  
##  Jun 12 - 22   : 3163  
##  Jun 23 - Jul 3: 2552  
##  Jul 4 - 13    : 2764

I also need to remove outliers, which are pretty standard erroneous points when using temperature-sensitive telemetry with devices external to the animal. I will follow Nicole Gaudenti’s 2021 paper, and exclude “any outliers greater than two standard deviations away from each lizard’s mean Tb”

Look at data distribution:

boxplot(temp_sub$ppTempC)

hist(temp_sub$ppTempC)

ggplot(temp_sub) +
  geom_point(aes(x = Serial,
                 y = ppTempC,
                 color = Serial)) +
  theme_classic()

ggplot(temp_sub) +
  geom_boxplot(aes(x = Serial,
                 y = ppTempC,
                 color = Serial)) +
  theme_classic()

Get mean and SD for each lizard’s Tb data, then use it to filter out extreme values:

temp_clean <- temp_sub %>%
  group_by(Serial) %>%
  #summarise(outs = boxplot.stats(ppTempC)$out) # less outliers found this way
  mutate(mean_Tb = mean(ppTempC),
         sd_Tb = sd(ppTempC)) %>%
  mutate(accept_low = mean_Tb - (2*sd_Tb),
         accept_high = mean_Tb + (2*sd_Tb)) %>%
  dplyr::filter(ppTempC < accept_high & ppTempC > accept_low) %>%
  dplyr::select(-mean_Tb, -sd_Tb, -accept_low, -accept_high)

summary(temp_clean)
##     DateTime                          Serial           Freq      
##  Min.   :2021-04-26 00:01:18.00   245-756: 4984   Min.   :150.0  
##  1st Qu.:2021-05-02 15:13:54.25   245-753: 4974   1st Qu.:150.4  
##  Median :2021-05-12 08:52:29.50   229-078: 4686   Median :151.1  
##  Mean   :2021-05-15 18:45:28.49   229-069: 4294   Mean   :151.0  
##  3rd Qu.:2021-05-22 09:58:02.75   252-886: 4251   3rd Qu.:151.5  
##  Max.   :2021-07-13 13:11:42.00   252-881: 3996   Max.   :152.0  
##                                   (Other):65645                  
##     ppTempC            pp         date_only           hour_of_day   
##  Min.   : 1.67   Min.   :1109   Min.   :2021-04-26   Min.   : 0.00  
##  1st Qu.:32.45   1st Qu.:1564   1st Qu.:2021-05-02   1st Qu.:10.00  
##  Median :36.48   Median :1694   Median :2021-05-12   Median :13.00  
##  Mean   :35.55   Mean   :1721   Mean   :2021-05-15   Mean   :12.75  
##  3rd Qu.:39.09   3rd Qu.:1862   3rd Qu.:2021-05-22   3rd Qu.:16.00  
##  Max.   :53.17   Max.   :2798   Max.   :2021-07-13   Max.   :23.00  
##                                                                     
##          date_chunks   
##  Apr 26 - May 6:36977  
##  May 10 - 20   :32519  
##  May 21 - 31   : 9080  
##  Jun 1 - 11    : 6005  
##  Jun 12 - 22   : 3067  
##  Jun 23 - Jul 3: 2494  
##  Jul 4 - 13    : 2688

What proportion of points did we remove?

(nrow(temp_sub) - nrow(temp_clean))/nrow(temp_sub)
## [1] 0.04548913

4.5%

Plot again:

ggplot(temp_clean) +
  geom_point(aes(x = Serial,
                 y = ppTempC,
                 color = Serial)) +
  theme_classic()

ggplot(temp_clean) +
  geom_boxplot(aes(x = Serial,
                 y = ppTempC,
                 color = Serial)) +
  theme_classic()

Individual IQRs definitely look reasonable.

Also, save the link between radio frequency and transmitter serial number to help with microhabitat data wrangling.

freq_serial <- temp_clean %>%
  group_by(Freq, Serial) %>%
  summarise(n = n()) #%>%
## `summarise()` has grouped output by 'Freq'. You can override using the
## `.groups` argument.
# check that there is 1 serial per frequency
#  group_by(Freq) %>%
 # summarise(n = n())

Microhabitat Use

This data was collected by radio-tracking lizards, finding them, then recording what type of microhabitat they were in. We have to read-in a few files, because at first the formatting was not consistent, then put the dataframes all together.

variables in this dataframe:

  • Serial = the radio-collar serial number for that observation (to be used to link to specific individuals)
  • DateTime = date and time the temperature was logged
  • microhabitat (was Micro_1) = which microhabitat lizards were observed in
  • Breed_Col = breeding coloration of lizards (we will not use)
  • bpmTempC = the temperature in Celsius from the lizards’ collar, based on the handheld receivers used to locate them, which record the pulse interval in beats per minute rather than the pulse period of the tower receiver for other Tb data
  • dwsTempOut = the temperature recorded by a weather station at the study site; these values were interpolated to the specific time of the observation since the weather station was only recording data every two hours ModTempC = the ground surface temperature estimated by a temperature model (Ian Axsom) for the location and time that each observation occurred; this is an estimate of the operative temperature available to these lizards

This is the nice, clean data that was collected early May onwards:

microhab_clean <- read.csv("./data/telemetry/behavior_observations.csv") %>%
  # force NAs on, then remove, observations for the "uncollared" lizard
  mutate(Serial_no = as.numeric(substr(Serial, 1, 5))) %>%
  dplyr::filter(complete.cases(Serial_no)) %>%
  # properly format data classes
  mutate(date_only = as.Date(substr(DateTime, 1, 10), format = "%Y-%m-%d"),
          hour_of_day = as.numeric(substr(DateTime, 12, 13)),
          DateTime = as.POSIXct(DateTime, format = "%Y-%m-%d %H:%M:%S"),
          Serial = factor(paste(paste(substr(as.character(Serial), 1, 3), 
                                     substr(as.character(Serial), 4, 7), sep = "-"))),
         microhabitat = factor(Micro_1,
                               labels = c("Burrow", # keep
                                          "Open", # change from Burrow slide
                                          "Open", # change from Burrow-Partial
                                          "Burrow", # change from Burrow/ephedra
                                          "Burrow", # change from Burrow/shrub
                                          "", ### change from Dead/Missing ###
                                          "Open", # change from Grass trims
                                          "", # change from IDK
                                          "", # change from No Visual
                                          "Open", # keep
                                          "Shade - Full", # keep
                                          "Shade - Partial", # keep
                                          "Shade - Partial", # remove notes
                                          "Shade - Partial", # remove notes
                                          "Shade - Partial", # remove notes
                                          "Shade - Partial", # remove notes
                                          "Shade - Partial", # remove notes
                                          "Shade - Partial", # remove notes
                                          "Shade - Partial" # remove notes
                                          ))
         ) %>%
  # for some reason, it only works to fix the order separately
  mutate(microhabitat = factor(microhabitat,
                               levels = c("Open", 
                                          "Shade - Partial", 
                                          "Shade - Full", 
                                          "Burrow"
                                          #"No Visual", "Dead/Missing"
                                          ))) %>%
  dplyr::select(Serial, DateTime, date_only, hour_of_day, microhabitat, 
                bpmTempC, dwsTempOut, ModTempC) %>%
  # add date interval chunk grouping notes
  dplyr::filter(date_only %nin% cap_dates)
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

This is the messy data collected in late April and very early May. Unfortunately, the following radio-frequencies were unable to be linked to a radio collar serial number: 150.14, 150.19, 150.21, and 150.28. But, that only applied to 11 observations.

microhab_early0 <- read.csv("./data/telemetry/behavior_observations_early_0.csv") %>%
  dplyr::select(date_only = date, time, 
                microhabitat = micro, 
                Serial = serial, 
                Freq = freq, freq_notes = X)
# this one must be split by whether freq or serial number of radio collars recorded
microhab_early0_freq_only <- microhab_early0 %>%
  dplyr::filter(!complete.cases(Serial)) %>%
  dplyr::select(-Serial)
microhab_early0_serial <- microhab_early0 %>%
  dplyr::filter(complete.cases(Serial)) %>%
  mutate(Serial = factor(paste(substr(Serial, 1, 3), 
                               substr(Serial, 4, 7), sep = "-")))

# load other csv, and add the dfs above to it
microhab_early_all <- read.csv("./data/telemetry/behavior_observations_early_1.csv") %>%
  dplyr::select(date_only = date, time, 
                microhabitat = micro, 
                Freq = freq, freq_notes) %>%
  rbind(microhab_early0_freq_only) %>%
  mutate(Freq = round(Freq, 2)) %>%
  left_join(freq_serial, by = 'Freq') %>%
  dplyr::select(-n) %>%
  rbind(microhab_early0_serial) %>%
  mutate(DateTime = as.POSIXct(paste(date_only, time, sep = " "), 
                               format = "%m/%d/%y %H:%M"),
         hour_of_day = as.numeric(substr(as.character(DateTime), 12, 13)),
         date_only = as.Date(date_only, format = "%m/%d/%y"),
         microhabitat = factor(str_to_title(microhabitat))) %>%
  dplyr::filter(complete.cases(Serial, microhabitat)) %>%
  # some of these were "focal observations" taken every few minutes
  # to be consistent with other obs, take first point of each focal
  group_by(Serial, date_only) %>%
  #summarise(min(date), max(date), min(time), max(time)) # check that data is one chunk
  # might remove a few usable data, but OK
  slice(1) %>%
  # remove these to make it easy to bind to clean one
  dplyr::select(-time, -Freq, -freq_notes) %>%
  # add these to make it easy to bind to clean one
  mutate(bpmTempC = NA, dwsTempOut = NA, ModTempC = NA)
summary(microhab_early_all) 
##    date_only                   microhabitat     Serial  
##  Min.   :2021-04-27   Burrow         :28    245-752: 5  
##  1st Qu.:2021-04-28   Open           :79    245-760: 5  
##  Median :2021-04-29   Shade - Full   : 4    252-884: 5  
##  Mean   :2021-04-29   Shade - Partial: 7    229-060: 4  
##  3rd Qu.:2021-05-01                         229-063: 4  
##  Max.   :2021-05-02                         229-069: 4  
##                                             (Other):91  
##     DateTime                       hour_of_day    bpmTempC       dwsTempOut    
##  Min.   :2021-05-01 10:03:00.00   Min.   : 2.00   Mode:logical   Mode:logical  
##  1st Qu.:2021-05-01 13:35:30.00   1st Qu.:10.50   NA's:118       NA's:118      
##  Median :2021-05-01 16:00:00.00   Median :13.00                                
##  Mean   :2021-05-01 22:00:40.35   Mean   :11.93                                
##  3rd Qu.:2021-05-02 11:00:30.00   3rd Qu.:15.00                                
##  Max.   :2021-05-02 15:46:00.00   Max.   :17.00                                
##  NA's   :63                       NA's   :63                                   
##  ModTempC      
##  Mode:logical  
##  NA's:118      
##                
##                
##                
##                
## 

Now, put together the original and messy-now-clean dataframes, and add some formatting:

microhab_pre <- microhab_clean %>%
  rbind(microhab_early_all) %>%
  dplyr::filter(complete.cases(microhabitat)) %>%
  left_join(group_dates, by = 'date_only') %>%
  left_join(liz_IDs_only, by = c('Serial' = 'radio_collar_serial')) %>%
  mutate(above_below = factor(case_when(microhabitat == "Burrow" ~ "Belowground",
                                        microhabitat != "Burrow" ~ "Aboveground"),
                              levels = c("Aboveground", "Belowground")))

summary(microhab_pre)
##       Serial        DateTime                        date_only         
##  245-744 : 100   Min.   :2021-05-01 10:03:00.00   Min.   :2021-04-27  
##  245-741 :  96   1st Qu.:2021-05-31 13:30:00.00   1st Qu.:2021-05-29  
##  229-078 :  93   Median :2021-06-22 13:17:00.00   Median :2021-06-21  
##  245-762a:  92   Mean   :2021-06-15 04:09:31.63   Mean   :2021-06-12  
##  229-079 :  91   3rd Qu.:2021-06-30 14:17:00.00   3rd Qu.:2021-06-30  
##  245-745 :  91   Max.   :2021-07-13 13:27:00.00   Max.   :2021-07-13  
##  (Other) :1213   NA's   :63                                           
##   hour_of_day             microhabitat    bpmTempC        dwsTempOut   
##  Min.   : 2.00   Open           :544   Min.   : 8.028   Min.   :10.11  
##  1st Qu.:10.00   Shade - Partial:147   1st Qu.:32.828   1st Qu.:26.61  
##  Median :12.00   Shade - Full   :214   Median :36.267   Median :31.28  
##  Mean   :11.99   Burrow         :871   Mean   :35.633   Mean   :30.25  
##  3rd Qu.:14.00                         3rd Qu.:39.792   3rd Qu.:34.60  
##  Max.   :18.00                         Max.   :49.823   Max.   :40.58  
##  NA's   :63                            NA's   :150      NA's   :173    
##     ModTempC             date_chunks  individual_ID    sex_M_F    
##  Min.   :21.24   Apr 26 - May 6:166   M-09   : 100   Female: 539  
##  1st Qu.:39.50   May 10 - 20   :158   F-10   :  96   Male  :1237  
##  Median :44.48   May 21 - 31   :181   M-14   :  93                
##  Mean   :43.77   Jun 1 - 11    :181   M-03a  :  92                
##  3rd Qu.:48.99   Jun 12 - 22   :253   M-08   :  91                
##  Max.   :57.86   Jun 23 - Jul 3:579   M-11   :  91                
##  NA's   :118     Jul 4 - 13    :258   (Other):1213                
##       above_below 
##  Aboveground:905  
##  Belowground:871  
##                   
##                   
##                   
##                   
## 

Check whether we have observations on lizards that were not alive or were not actually wearing their radio-collars because they dropped them.

microhab_check <- microhab_pre %>%
  left_join(check_gone, by = 'Serial') %>%
  dplyr::filter(date_only > date_def_gone) %>%
  arrange(Serial) %>%
  dplyr::select(individual_ID, Serial, date_only, date_def_gone)
microhab_check
##   individual_ID  Serial  date_only date_def_gone
## 1          F-17 229-060 2021-05-19    2021-05-18
## 2          F-08 229-070 2021-05-11    2021-05-07
## 3          F-08 229-070 2021-05-23    2021-05-07
## 4          F-08 229-070 2021-06-07    2021-05-07
## 5          M-03 245-762 2021-06-10    2021-05-07

Check what the name switch should be for F-08 and M-03:

unique(liz_dat$radio_collar_serial)
##  [1] 245-761  252-886  245-762  252-887  245-753  245-752  245-754  229-079 
##  [9] 245-744  229-069  252-885  245-759  245-751  245-757  229-059  229-063 
## [17] 245-747  229-070  229-065  245-741  245-748  252-882  252-881  229-078 
## [25] 245-750  229-080  245-760  245-746  245-745  245-742  252-884  252-883 
## [33] 245-756  229-073  229-066  229-072  229-060  229-074  229-077  245-762a
## [41] 229-070a
## 41 Levels: 229-059 229-060 229-063 229-065 229-066 229-069 229-070 ... 252-887
unique(liz_dat$individual_ID)
##  [1] M-01  M-02  M-03  M-04  M-05  M-06  M-07  M-08  M-09  M-10  F-01  F-02 
## [13] F-03  F-04  F-05  F-06  F-07  F-08  F-09  F-10  M-11  M-12  M-13  M-14 
## [25] M-15  M-16  M-17  M-18  M-19  M-20  F-11  F-12  F-13  F-14  F-15  F-16 
## [37] F-17  F-18  F-19  M-03a F-08a
## 79 Levels: F-01 F-02 F-03 F-04 F-05 F-06 F-07 F-08 F-08a F-09 F-10 ... W-039

Then fix those things:

microhab <- microhab_pre %>%
  # delete the one obs for F-17 that was outside the valid date range
  dplyr::filter(!(individual_ID == "F-17" & date_only == "2021-05-19")) %>%
  # change the serial and ID of F-08 to be "-a" bc collar was switched
          # for F-08a
  mutate(individual_ID = replace(individual_ID, # what I want to change
                                 individual_ID == "F-08" & # if this
                                   date_only > "2021-05-07", # and this
                                 "F-08a"), # change to this
         Serial = replace(Serial, # what I want to change
                          Serial == "229-070" & # if this
                          date_only > "2021-05-07", # and this
                          "229-070a"), # change to this
         # for M-03a
         individual_ID = replace(individual_ID, # what I want to change
                                 individual_ID == "M-03" & # if this
                                   date_only > "2021-05-07", # and this
                                 "M-03a"), # change to this
         Serial = replace(Serial, # what I want to change
                          Serial == "245-762" & # if this
                          date_only > "2021-05-07", # and this
                          "245-762a") # change to this
         )

Check that worked:

nrow(microhab) == nrow(microhab_pre) -1
## [1] TRUE
microhab_check2 <- microhab %>%
  left_join(check_gone, by = 'Serial') %>%
  dplyr::filter(date_only > date_def_gone) %>%
  arrange(Serial) %>%
  dplyr::select(individual_ID, Serial, date_only, date_def_gone)
microhab_check2
## [1] individual_ID Serial        date_only     date_def_gone
## <0 rows> (or 0-length row.names)

Export full microhab df

#write_rds(microhab, "./data/telemetry_clean_2021.RDS")

methods stats: How many observations and days was each lizard tracked for?

micro_obs_ns <- microhab %>%
  group_by(Serial) %>%
  summarise(n = n(),
            min_date = min(date_only),
            max_date = max(date_only),
            date_range = max_date - min_date) %>%
  mutate(per_day = n/as.numeric(date_range),
         per_wk = per_day*7)
micro_obs_ns %>%
  dplyr::filter(Serial != "252-882") %>%
  summarise(mean(n), sd(n),
            mean(date_range), sd(date_range),
            mean(per_wk), sd(per_wk))
## # A tibble: 1 × 6
##   `mean(n)` `sd(n)` `mean(date_range)` `sd(date_range)` `mean(per_wk)` sd(per_…¹
##       <dbl>   <dbl> <drtn>                        <dbl>          <dbl>     <dbl>
## 1      46.7    38.5 47.13158 days                  28.8           6.09      2.50
## # … with abbreviated variable name ¹​`sd(per_wk)`

Thermal Metrics

For each lizard, we want to calculate some assessment of how hot they let themselves be, and their thermoregulatory accuracy.

Daily Maximums

daily_max <- temp_clean %>%
  # quick look to use
  #dplyr::filter(date_chunks %in% c("Apr 26 - May 6", "May 10 - 20")) %>%
  
  # first do daily maximums
  group_by(date_only, date_chunks, Serial) %>%
  mutate(Tb_percentile_50 = quantile(ppTempC, 0.5),
         Tb_percentile_60 = quantile(ppTempC, 0.6),
         Tb_percentile_70 = quantile(ppTempC, 0.7),
         Tb_percentile_80 = quantile(ppTempC, 0.8),
         Tb_percentile_90 = quantile(ppTempC, 0.9)) %>%
  ungroup() %>%
  # select only the data between the 80-90th percentiles of Tb for each lizard
  dplyr::filter(ppTempC > Tb_percentile_80 & ppTempC < Tb_percentile_90) %>%
  # calculate the mean Tb for that subset for each lizard per day
  group_by(date_only, date_chunks, Serial) %>%
  summarise(daily_Tb_percentile_80_90 = mean(ppTempC)) %>%
  # then do the mean Tb between 80-90%ile for each lizard per time chunk (~11 days)
  group_by(date_chunks, Serial) %>%
  summarise(chunk_avg_Tb_percentile_80_90 = mean(daily_Tb_percentile_80_90))
## `summarise()` has grouped output by 'date_only', 'date_chunks'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.
summary(daily_max)
##          date_chunks     Serial    chunk_avg_Tb_percentile_80_90
##  Apr 26 - May 6:35   229-065:  7   Min.   :18.89                
##  May 10 - 20   :32   229-069:  7   1st Qu.:37.51                
##  May 21 - 31   :28   229-078:  7   Median :39.38                
##  Jun 1 - 11    :23   245-742:  7   Mean   :38.37                
##  Jun 12 - 22   :18   245-748:  7   3rd Qu.:40.79                
##  Jun 23 - Jul 3:15   245-753:  7   Max.   :46.85                
##  Jul 4 - 13    :12   (Other):121

Db

From Kat Ivey’s 2020 paper, Tset range is: 32.3±1.2 - 37.5±1.1 deg C.

From Nicole Gaudenti’s 2021 paper, Tset range: 33.2-37.9 deg C.

I will use the widest margin created by both: 32.3 (Ivey) - 37.9 (Gaudenti)

Db <- temp_clean %>%
  # we only want daytime Tb for this
  # currently using Kat's limits
  dplyr::filter(hour_of_day <=19 & hour_of_day >=7) %>%
  mutate(Tpref_low = 32.3,
         Tpref_high = 37.9) %>%
  # calculate Db for every Tb measurement
  mutate(Db = case_when(#ppTempC < Tpref_low ~ ppTempC - Tpref_low, # for pos/neg version
                        ppTempC < Tpref_low ~ Tpref_low - ppTempC, # for abs value version
                        ppTempC > Tpref_high ~ ppTempC - Tpref_high,
                        ppTempC < Tpref_high & ppTempC > Tpref_low ~ 0)) %>%
  # then avg Db by lizard & date chunk
  group_by(date_chunks, Serial) %>%
  summarise(chunk_mean_Db_deg_C_diff = mean(Db))
## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.
summary(Db)
##          date_chunks     Serial    chunk_mean_Db_deg_C_diff
##  Apr 26 - May 6:35   229-065:  7   Min.   : 0.000          
##  May 10 - 20   :32   229-069:  7   1st Qu.: 1.033          
##  May 21 - 31   :29   229-078:  7   Median : 1.532          
##  Jun 1 - 11    :25   229-079:  7   Mean   : 2.453          
##  Jun 12 - 22   :20   245-742:  7   3rd Qu.: 2.303          
##  Jun 23 - Jul 3:18   245-745:  7   Max.   :16.429          
##  Jul 4 - 13    :16   (Other):133

Microhabitat Use

Proportions of Observations

To calculate microhabitat use, I literally just calculate the proportion of observations for a given individual/time range that were in each microhabitat.

So, for each time chunk, for each lizard (radio-collar serial), I need to get the following values:

  • number of total observations
  • number of observations in each microhabitat
  • proportion of observations in each microhabitat
  • aboveground temps based on Ian’s model ?

We do not have a lot of obs for each lizard for every time chunk, so how many observations should there be to make the proportions valid/reasonable? One potential way to do this could be to at least have 1-2 obs for every microhabitat type we are trying to observe (so, minimum points per date chunk either 4 or 8, once I get rid of the irrelevant ones).

# where to cut off n obs?
microhab %>%
  # get only 1 per lizard
  group_by(date_chunks, Serial) %>%
  summarise(n_obs = n()) %>%
  ungroup() %>%
  # How many lizards have each (n) obs?
  group_by(n_obs) %>%
  summarise(n = n())
## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.
## # A tibble: 26 × 2
##    n_obs     n
##    <int> <int>
##  1     1     5
##  2     2    12
##  3     3     6
##  4     4    18
##  5     5    27
##  6     6    18
##  7     7    13
##  8     8     4
##  9     9     7
## 10    10     9
## # … with 16 more rows

It seems like a good cut-off is a minimum of 4 observations per lizard, per time chunk. Otherwise, we would lose a LOT of data.

# micro use BY individual lizard
MH_use_by_indiv <- microhab %>%
  group_by(Serial, date_chunks, microhabitat) %>%
  summarise(n_obs_micro = n()) %>%
  group_by(Serial, date_chunks) %>%
  mutate(n_obs_total = sum(n_obs_micro)) %>%
  ungroup() %>%
  # expect that this takes out 42 obs/rows
  dplyr::filter(n_obs_total >= 4) %>%
  mutate(MH_use_proportion = n_obs_micro/n_obs_total) #%>%
## `summarise()` has grouped output by 'Serial', 'date_chunks'. You can override
## using the `.groups` argument.
# check that all proportions for a given lizard and time chunk = 1
#  group_by(Serial, date_chunks) %>%
 # summarise(prop_total = sum(MH_use_proportion))
head(MH_use_by_indiv)
## # A tibble: 6 × 6
##   Serial  date_chunks    microhabitat    n_obs_micro n_obs_total MH_use_propor…¹
##   <fct>   <fct>          <fct>                 <int>       <int>           <dbl>
## 1 229-059 Apr 26 - May 6 Open                      3           4           0.75 
## 2 229-059 Apr 26 - May 6 Burrow                    1           4           0.25 
## 3 229-060 Apr 26 - May 6 Open                      4           6           0.667
## 4 229-060 Apr 26 - May 6 Burrow                    2           6           0.333
## 5 229-060 May 10 - 20    Open                      2           4           0.5  
## 6 229-060 May 10 - 20    Shade - Partial           1           4           0.25 
## # … with abbreviated variable name ¹​MH_use_proportion
summary(MH_use_by_indiv)
##       Serial            date_chunks          microhabitat  n_obs_micro    
##  245-741 : 23   Apr 26 - May 6:62   Open           :122   Min.   : 1.000  
##  229-079 : 22   May 10 - 20   :67   Shade - Partial: 75   1st Qu.: 1.000  
##  245-745 : 22   May 21 - 31   :70   Shade - Full   : 67   Median : 3.000  
##  245-744 : 21   Jun 1 - 11    :59   Burrow         :118   Mean   : 4.524  
##  245-762a: 19   Jun 12 - 22   :50                         3rd Qu.: 5.000  
##  245-761 : 18   Jun 23 - Jul 3:45                         Max.   :32.000  
##  (Other) :257   Jul 4 - 13    :29                                         
##   n_obs_total    MH_use_proportion
##  Min.   : 4.00   Min.   :0.0303   
##  1st Qu.: 5.00   1st Qu.:0.2000   
##  Median : 8.00   Median :0.3077   
##  Mean   :11.56   Mean   :0.4005   
##  3rd Qu.:14.00   3rd Qu.:0.6000   
##  Max.   :37.00   Max.   :1.0000   
## 
# micro use for all lizards in the study, aggregated
MH_use_across_indivs <- microhab %>%
  group_by(date_chunks, microhabitat) %>%
  summarise(n_obs_micro = n(),
            mean_mod_temp = mean(ModTempC, na.rm = T)) %>%
  group_by(date_chunks) %>%
  mutate(n_obs_total = sum(n_obs_micro),
         mean_mod_temp2 = mean(mean_mod_temp, na.rm = T)) %>%
  ungroup() %>%
  mutate(MH_use_proportion = n_obs_micro/n_obs_total)
## `summarise()` has grouped output by 'date_chunks'. You can override using the
## `.groups` argument.
head(MH_use_across_indivs)
## # A tibble: 6 × 7
##   date_chunks    microhabitat    n_obs_micro mean_mod_…¹ n_obs…² mean_…³ MH_us…⁴
##   <fct>          <fct>                 <int>       <dbl>   <int>   <dbl>   <dbl>
## 1 Apr 26 - May 6 Open                    110        37.3     166    40.5  0.663 
## 2 Apr 26 - May 6 Shade - Partial          13        42.7     166    40.5  0.0783
## 3 Apr 26 - May 6 Shade - Full              8        41.4     166    40.5  0.0482
## 4 Apr 26 - May 6 Burrow                   35        40.5     166    40.5  0.211 
## 5 May 10 - 20    Open                     82        42.3     157    44.3  0.522 
## 6 May 10 - 20    Shade - Partial          23        44.3     157    44.3  0.146 
## # … with abbreviated variable names ¹​mean_mod_temp, ²​n_obs_total,
## #   ³​mean_mod_temp2, ⁴​MH_use_proportion
# micro use for all lizards, BY SEX
MH_use_across_indivs_sex <- microhab %>%
  group_by(date_chunks, microhabitat, sex_M_F) %>%
  summarise(n_obs_micro = n(),
            mean_mod_temp = mean(ModTempC, na.rm = T)) %>%
  group_by(date_chunks, sex_M_F) %>%
  mutate(n_obs_total = sum(n_obs_micro),
         mean_mod_temp2 = mean(mean_mod_temp, na.rm = T)) %>%
  ungroup() %>%
  mutate(MH_use_proportion = n_obs_micro/n_obs_total)
## `summarise()` has grouped output by 'date_chunks', 'microhabitat'. You can
## override using the `.groups` argument.
summary(MH_use_across_indivs_sex)
##          date_chunks          microhabitat   sex_M_F    n_obs_micro    
##  Apr 26 - May 6:8    Open           :14    Female:28   Min.   :  1.00  
##  May 10 - 20   :8    Shade - Partial:14    Male  :28   1st Qu.:  9.00  
##  May 21 - 31   :8    Shade - Full   :14                Median : 17.00  
##  Jun 1 - 11    :8    Burrow         :14                Mean   : 31.70  
##  Jun 12 - 22   :8                                      3rd Qu.: 31.75  
##  Jun 23 - Jul 3:8                                      Max.   :269.00  
##  Jul 4 - 13    :8                                                      
##  mean_mod_temp    n_obs_total    mean_mod_temp2  MH_use_proportion
##  Min.   :33.22   Min.   : 52.0   Min.   :39.76   Min.   :0.01724  
##  1st Qu.:40.21   1st Qu.: 69.0   1st Qu.:40.56   1st Qu.:0.08974  
##  Median :43.50   Median : 86.5   Median :43.75   Median :0.15281  
##  Mean   :43.08   Mean   :126.8   Mean   :43.08   Mean   :0.25000  
##  3rd Qu.:45.85   3rd Qu.:141.0   3rd Qu.:44.92   3rd Qu.:0.37911  
##  Max.   :49.76   Max.   :438.0   Max.   :45.23   Max.   :0.77000  
## 

Proportion of time aboveground

A version of microhabitat use that’s basically burrow/not.

prop_above_by_indiv <- MH_use_by_indiv %>%
  dplyr::filter(microhabitat == "Burrow") %>%
  mutate(prop_above = 1 - MH_use_proportion)
prop_above_across_indiv <- MH_use_across_indivs %>%
  dplyr::filter(microhabitat == "Burrow") %>%
  mutate(prop_above = 1 - MH_use_proportion)
prop_above_across_indiv_by_sex <- MH_use_across_indivs_sex %>%
  dplyr::filter(microhabitat == "Burrow") %>%
  mutate(prop_above = 1 - MH_use_proportion)

Onset of Estivation

This didn’t make it into the paper, because the estivation onset dates we got were somewhat questionable, and we have no way of knowing whether it is due to our tracking methods and the dropping/loss of collars throughout the season.

# first, get some key dates for each lizard
last_date <- microhab %>%
  # for each lizard
  group_by(Serial) %>%
  # what was the latest date we tracked them?
  summarise(date_last_obs = max(date_only),
            # and how many observations did they have total?
            n_obs_total = n())

last_date_above <- microhab %>%
  # first, subset aboveground obs only
  dplyr::filter(microhabitat != "Burrow") %>%
  # for each lizard
  group_by(Serial) %>%
  # what was the latest date we observed each lizard ABOVEground?
  summarise(est_onset = max(date_only))
  
estivating <- microhab %>%
  left_join(last_date_above, by = 'Serial') %>%
  dplyr::filter(date_only > est_onset) %>%
  group_by(Serial, est_onset) %>%
  # how many times did we observe each lizard while they were estivating?
  summarise(n_obs_est = n()) %>%
  left_join(last_date, by = 'Serial') %>%
  mutate(est_obs_days = as.numeric(date_last_obs - est_onset),
         prop_obs_est = n_obs_est/n_obs_total) %>%
  # make sure we tracked them after "last" aboveground observation
  dplyr::filter(est_onset < date_last_obs) %>%
  # and we still had a signal on them till the end
  dplyr::filter(date_last_obs == "2021-07-13")
## `summarise()` has grouped output by 'Serial'. You can override using the
## `.groups` argument.
summary(estivating)
##       Serial    est_onset            n_obs_est     date_last_obs       
##  229-063 :1   Min.   :2021-05-11   Min.   : 7.00   Min.   :2021-07-13  
##  229-069 :1   1st Qu.:2021-06-11   1st Qu.:30.25   1st Qu.:2021-07-13  
##  229-070a:1   Median :2021-06-20   Median :48.00   Median :2021-07-13  
##  229-078 :1   Mean   :2021-06-18   Mean   :42.29   Mean   :2021-07-13  
##  229-079 :1   3rd Qu.:2021-06-27   3rd Qu.:58.75   3rd Qu.:2021-07-13  
##  245-742 :1   Max.   :2021-07-08   Max.   :65.00   Max.   :2021-07-13  
##  (Other) :8                                                            
##   n_obs_total     est_obs_days    prop_obs_est    
##  Min.   :70.00   Min.   : 5.00   Min.   :0.07692  
##  1st Qu.:83.50   1st Qu.:15.75   1st Qu.:0.38155  
##  Median :86.50   Median :22.50   Median :0.55868  
##  Mean   :85.57   Mean   :24.43   Mean   :0.50250  
##  3rd Qu.:90.75   3rd Qu.:31.25   3rd Qu.:0.65830  
##  Max.   :93.00   Max.   :63.00   Max.   :0.92857  
## 
estivating %>%
  mutate(est = "Yes") %>%
  group_by(est) %>%
  summarise(mean(n_obs_est),
            sd(n_obs_est),
            min(n_obs_est),
            max(n_obs_est))
## # A tibble: 1 × 5
##   est   `mean(n_obs_est)` `sd(n_obs_est)` `min(n_obs_est)` `max(n_obs_est)`
##   <chr>             <dbl>           <dbl>            <int>            <int>
## 1 Yes                42.3            19.8                7               65

Max Temps ~ Hydric

Now, look for relationships between hydric physiology and thermoregulation.

Prep Data

max_temps_hydric <- liz_dat_apr_may_only %>%
  left_join(daily_max, by = c('date_chunks', 
                              'radio_collar_serial' = 'Serial')) %>%
  # remove NAs
  dplyr::filter(complete.cases(chunk_avg_Tb_percentile_80_90)) %>%
  # remove the one outlier wayyyyy below the other values
  dplyr::filter(chunk_avg_Tb_percentile_80_90>30)

summary(max_temps_hydric)
##  capture_date_time              capture_date        radio_collar_serial
##  Min.   :2021-04-23 03:39:00   Min.   :2021-04-23   229-059: 2         
##  1st Qu.:2021-04-23 07:06:30   1st Qu.:2021-04-23   229-060: 2         
##  Median :2021-04-24 03:49:00   Median :2021-04-24   229-063: 2         
##  Mean   :2021-04-28 03:39:18   Mean   :2021-04-28   229-069: 2         
##  3rd Qu.:2021-05-07 04:27:00   3rd Qu.:2021-05-07   229-072: 2         
##  Max.   :2021-05-08 07:19:00   Max.   :2021-05-08   229-077: 2         
##  NA's   :4                                          (Other):42         
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:54          F-01   : 2    Female:27   Min.   :26.00   Min.   : 85.0  
##  Class :character   F-05   : 2    Male  :27   1st Qu.:33.17   1st Qu.: 95.0  
##  Mode  :character   F-06   : 2                Median :37.50   Median : 99.0  
##                     F-11   : 2                Mean   :38.93   Mean   :100.4  
##                     F-12   : 2                3rd Qu.:43.75   3rd Qu.:107.8  
##                     F-13   : 2                Max.   :53.60   Max.   :122.0  
##                     (Other):42                                               
##         tmt     tmt_mass_change_g capture_to_msmt   hematocrit_percent
##  Water Tmt:25   Min.   :-3.0000   Min.   :  11.00   Min.   :23.00     
##  Sham Tmt :29   1st Qu.:-1.0000   1st Qu.:  59.75   1st Qu.:30.00     
##  No Tmt   : 0   Median : 0.0000   Median : 152.00   Median :36.00     
##                 Mean   :-0.3111   Mean   : 410.56   Mean   :34.81     
##                 3rd Qu.: 0.0000   3rd Qu.: 265.00   3rd Qu.:39.00     
##                 Max.   : 3.0000   Max.   :1665.00   Max.   :58.00     
##                                   NA's   :4         NA's   :1         
##  cloacal_temp_C   CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :24.0   Min.   : 0.650   Min.   :19.07   Min.   :12.20  
##  1st Qu.:32.0   1st Qu.: 7.077   1st Qu.:28.71   1st Qu.:14.53  
##  Median :33.0   Median : 9.262   Median :30.62   Median :17.87  
##  Mean   :32.5   Mean   : 9.481   Mean   :29.49   Mean   :19.30  
##  3rd Qu.:34.0   3rd Qu.:11.377   3rd Qu.:31.86   3rd Qu.:22.38  
##  Max.   :36.0   Max.   :21.673   Max.   :33.55   Max.   :35.83  
##                                                                 
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.205   Min.   :1.415   Min.   :316.0           April:34   1 :34  
##  1st Qu.:3.938   1st Qu.:3.064   1st Qu.:347.7           May  :20   3 :20  
##  Median :4.396   Median :3.694   Median :365.0           July : 0   12: 0  
##  Mean   :4.188   Mean   :3.421   Mean   :367.3                             
##  3rd Qu.:4.718   3rd Qu.:4.035   3rd Qu.:381.3                             
##  Max.   :5.188   Max.   :4.319   Max.   :436.5                             
##                                  NA's   :1                                 
##     week_num     ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1.000   Min.   :2021-04-24   Gravid    :14   Min.   :3.000  
##  1st Qu.:1.000   1st Qu.:2021-04-24   Not Gravid:11   1st Qu.:3.250  
##  Median :1.000   Median :2021-04-24   NA's      :29   Median :4.000  
##  Mean   :1.741   Mean   :2021-04-29                   Mean   :3.714  
##  3rd Qu.:3.000   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3.000   Max.   :2021-05-08                   Max.   :4.000  
##                  NA's   :29                           NA's   :40     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage   
##  Min.   :0.6400          Min.   :1.820                small round: 1  
##  1st Qu.:0.8775          1st Qu.:2.502                large round: 8  
##  Median :0.9950          Median :2.885                soft       : 2  
##  Mean   :1.1093          Mean   :3.140                soft oblong: 0  
##  3rd Qu.:1.3750          3rd Qu.:3.970                firm oblong: 2  
##  Max.   :1.8000          Max.   :4.740                hard oblong: 1  
##  NA's   :40              NA's   :40                   NA's       :40  
##    dev_point          SMI                date_chunks
##  Min.   :1.000   Min.   :27.04   Apr 26 - May 6:34  
##  1st Qu.:2.000   1st Qu.:33.00   May 10 - 20   :20  
##  Median :2.000   Median :34.73   May 21 - 31   : 0  
##  Mean   :2.714   Mean   :35.52   Jun 1 - 11    : 0  
##  3rd Qu.:3.000   3rd Qu.:38.44   Jun 12 - 22   : 0  
##  Max.   :5.000   Max.   :45.63   Jun 23 - Jul 3: 0  
##  NA's   :40                      Jul 4 - 13    : 0  
##  chunk_avg_Tb_percentile_80_90
##  Min.   :35.53                
##  1st Qu.:38.35                
##  Median :39.14                
##  Mean   :39.05                
##  3rd Qu.:39.81                
##  Max.   :41.64                
## 

Plots

sub_maxs <- max_temps_hydric %>%
  dplyr::select(chunk_avg_Tb_percentile_80_90, CEWL_g_m2h, osmolality_mmol_kg_mean, mass_g, SMI)
PerformanceAnalytics::chart.Correlation(sub_maxs)
## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

ggplot(max_temps_hydric) +
  geom_point(aes(x = CEWL_g_m2h,
                 y = chunk_avg_Tb_percentile_80_90,
                 color = date_chunks)) +
  geom_smooth(aes(x = CEWL_g_m2h,
                 y = chunk_avg_Tb_percentile_80_90),
              se = F,
              method = "lm",
              formula = y~x)

ggplot(max_temps_hydric) +
  geom_point(aes(x = osmolality_mmol_kg_mean,
                 y = chunk_avg_Tb_percentile_80_90,
                 color = date_chunks)) +
  geom_smooth(aes(x = osmolality_mmol_kg_mean,
                 y = chunk_avg_Tb_percentile_80_90),
              se = F,
              method = "lm",
              formula = y~x)
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

Models

# CEWL
max_temp_CEWL_mod <- lmerTest::lmer(data = max_temps_hydric,
                                chunk_avg_Tb_percentile_80_90 ~ CEWL_g_m2h +
                                  (1|individual_ID))
summary(max_temp_CEWL_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ CEWL_g_m2h + (1 | individual_ID)
##    Data: max_temps_hydric
## 
## REML criterion at convergence: 174.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.07592 -0.47459 -0.06368  0.51363  1.69421 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.6737   0.8208  
##  Residual                  0.7941   0.8911  
## Number of obs: 54, groups:  individual_ID, 36
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 39.13941    0.38430 48.59891 101.845   <2e-16 ***
## CEWL_g_m2h  -0.02071    0.03586 39.07870  -0.578    0.567    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## CEWL_g_m2h -0.876
#anova(max_temp_CEWL_mod, type = "3", ddf = "Kenward-Roger")
plot(max_temp_CEWL_mod)

# osmolality
max_temp_osml_mod <- lmerTest::lmer(data = max_temps_hydric,
                                chunk_avg_Tb_percentile_80_90 ~ osmolality_mmol_kg_mean +
                                  (1|individual_ID))
summary(max_temp_osml_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ osmolality_mmol_kg_mean + (1 |  
##     individual_ID)
##    Data: max_temps_hydric
## 
## REML criterion at convergence: 172.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0620 -0.5175 -0.0073  0.5198  1.6943 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.7480   0.8649  
##  Residual                  0.6931   0.8325  
## Number of obs: 53, groups:  individual_ID, 35
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)             35.262208   2.058746 45.768535  17.128   <2e-16 ***
## osmolality_mmol_kg_mean  0.010069   0.005586 45.564647   1.803   0.0781 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.996
#anova(max_temp_osml_mod, type = "3", ddf = "Kenward-Roger")
plot(max_temp_osml_mod)

# mass
max_temp_mass_mod <- lmerTest::lmer(data = max_temps_hydric,
                                chunk_avg_Tb_percentile_80_90 ~ mass_g +
                                  (1|individual_ID))
summary(max_temp_mass_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ mass_g + (1 | individual_ID)
##    Data: max_temps_hydric
## 
## REML criterion at convergence: 173.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.16856 -0.44986 -0.05865  0.47082  1.81093 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.6622   0.8138  
##  Residual                  0.7724   0.8789  
## Number of obs: 54, groups:  individual_ID, 36
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 37.83524    0.89481 40.51821  42.283   <2e-16 ***
## mass_g       0.02870    0.02266 40.04802   1.267    0.213    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## mass_g -0.979
#anova(max_temp_mass_mod, type = "3", ddf = "Kenward-Roger")
plot(max_temp_mass_mod)

# SMI
max_temp_SMI_mod <- lmerTest::lmer(data = max_temps_hydric,
                                chunk_avg_Tb_percentile_80_90 ~ SMI +
                                  (1|individual_ID))
summary(max_temp_SMI_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ SMI + (1 | individual_ID)
##    Data: max_temps_hydric
## 
## REML criterion at convergence: 169.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15402 -0.50939 -0.09287  0.40471  1.99940 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.5610   0.7490  
##  Residual                  0.7738   0.8797  
## Number of obs: 54, groups:  individual_ID, 36
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 36.18079    1.29543 40.90388   27.93   <2e-16 ***
## SMI          0.07826    0.03623 40.50454    2.16   0.0368 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## SMI -0.991
#anova(max_temp_SMI_mod, type = "3", ddf = "Kenward-Roger")
plot(max_temp_SMI_mod)

# supplemental hydration tmt
max_temp_tmt_mod <- lmerTest::lmer(data = max_temps_hydric,
                                chunk_avg_Tb_percentile_80_90 ~ tmt +
                                  (1|individual_ID))
summary(max_temp_tmt_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_avg_Tb_percentile_80_90 ~ tmt + (1 | individual_ID)
##    Data: max_temps_hydric
## 
## REML criterion at convergence: 166.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0046 -0.5040 -0.0538  0.5567  1.6822 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.5918   0.7693  
##  Residual                  0.7909   0.8893  
## Number of obs: 54, groups:  individual_ID, 36
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  38.6460     0.2542 35.9445 152.052   <2e-16 ***
## tmtSham Tmt   0.6050     0.3578 31.6304   1.691    0.101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tmtSham Tmt -0.710
#anova(max_temp_tmt_mod, type = "3", ddf = "Kenward-Roger")
plot(max_temp_tmt_mod)

Db ~ Hydric

Prep Data

db_hydric <- liz_dat_apr_may_only %>%
  left_join(Db, by = c('date_chunks', 
                              'radio_collar_serial' = 'Serial')) %>%
  # remove NAs
  dplyr::filter(complete.cases(chunk_mean_Db_deg_C_diff)) %>%
  # remove outliers
  dplyr::filter(chunk_mean_Db_deg_C_diff > -4)

summary(db_hydric)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 03:39:00.00   Min.   :2021-04-23   229-059: 2         
##  1st Qu.:2021-04-23 07:13:00.00   1st Qu.:2021-04-23   229-060: 2         
##  Median :2021-04-24 03:49:00.00   Median :2021-04-24   229-063: 2         
##  Mean   :2021-04-28 01:49:02.34   Mean   :2021-04-28   229-069: 2         
##  3rd Qu.:2021-05-07 04:23:00.00   3rd Qu.:2021-05-07   229-072: 2         
##  Max.   :2021-05-08 07:19:00.00   Max.   :2021-05-08   229-077: 2         
##  NA's   :4                                             (Other):43         
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:55          F-01   : 2    Female:28   Min.   :26.00   Min.   : 85.0  
##  Class :character   F-05   : 2    Male  :27   1st Qu.:33.35   1st Qu.: 95.0  
##  Mode  :character   F-06   : 2                Median :38.00   Median : 99.0  
##                     F-11   : 2                Mean   :38.91   Mean   :100.6  
##                     F-12   : 2                3rd Qu.:43.60   3rd Qu.:108.0  
##                     F-13   : 2                Max.   :53.60   Max.   :122.0  
##                     (Other):43                                               
##         tmt     tmt_mass_change_g capture_to_msmt  hematocrit_percent
##  Water Tmt:26   Min.   :-3.0000   Min.   :  11.0   Min.   :23.00     
##  Sham Tmt :29   1st Qu.:-1.0000   1st Qu.:  60.5   1st Qu.:30.00     
##  No Tmt   : 0   Median : 0.0000   Median : 150.0   Median :36.00     
##                 Mean   :-0.3236   Mean   : 403.8   Mean   :34.63     
##                 3rd Qu.: 0.0000   3rd Qu.: 264.0   3rd Qu.:38.75     
##                 Max.   : 3.0000   Max.   :1665.0   Max.   :58.00     
##                                   NA's   :4        NA's   :1         
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :24.00   Min.   : 0.650   Min.   :19.07   Min.   :12.20  
##  1st Qu.:32.00   1st Qu.: 7.377   1st Qu.:28.74   1st Qu.:14.53  
##  Median :33.00   Median : 9.318   Median :30.68   Median :18.67  
##  Mean   :32.53   Mean   : 9.649   Mean   :29.54   Mean   :19.32  
##  3rd Qu.:34.00   3rd Qu.:11.886   3rd Qu.:31.88   3rd Qu.:22.30  
##  Max.   :36.00   Max.   :21.673   Max.   :33.55   Max.   :35.83  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.205   Min.   :1.415   Min.   :316.0           April:35   1 :35  
##  1st Qu.:3.945   1st Qu.:3.068   1st Qu.:346.0           May  :20   3 :20  
##  Median :4.411   Median :3.703   Median :364.2           July : 0   12: 0  
##  Mean   :4.199   Mean   :3.428   Mean   :366.9                             
##  3rd Qu.:4.723   3rd Qu.:4.033   3rd Qu.:381.2                             
##  Max.   :5.188   Max.   :4.319   Max.   :436.5                             
##                                  NA's   :1                                 
##     week_num     ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1.000   Min.   :2021-04-24   Gravid    :14   Min.   :3.000  
##  1st Qu.:1.000   1st Qu.:2021-04-24   Not Gravid:12   1st Qu.:3.250  
##  Median :1.000   Median :2021-04-24   NA's      :29   Median :4.000  
##  Mean   :1.727   Mean   :2021-04-29                   Mean   :3.714  
##  3rd Qu.:3.000   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3.000   Max.   :2021-05-08                   Max.   :4.000  
##                  NA's   :29                           NA's   :41     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage   
##  Min.   :0.6400          Min.   :1.820                small round: 1  
##  1st Qu.:0.8775          1st Qu.:2.502                large round: 8  
##  Median :0.9950          Median :2.885                soft       : 2  
##  Mean   :1.1093          Mean   :3.140                soft oblong: 0  
##  3rd Qu.:1.3750          3rd Qu.:3.970                firm oblong: 2  
##  Max.   :1.8000          Max.   :4.740                hard oblong: 1  
##  NA's   :41              NA's   :41                   NA's       :41  
##    dev_point          SMI                date_chunks chunk_mean_Db_deg_C_diff
##  Min.   :1.000   Min.   :27.04   Apr 26 - May 6:35   Min.   : 0.5738         
##  1st Qu.:2.000   1st Qu.:33.00   May 10 - 20   :20   1st Qu.: 1.0201         
##  Median :2.000   Median :34.71   May 21 - 31   : 0   Median : 1.2748         
##  Mean   :2.714   Mean   :35.37   Jun 1 - 11    : 0   Mean   : 1.7030         
##  3rd Qu.:3.000   3rd Qu.:38.43   Jun 12 - 22   : 0   3rd Qu.: 1.6024         
##  Max.   :5.000   Max.   :45.63   Jun 23 - Jul 3: 0   Max.   :11.5644         
##  NA's   :41                      Jul 4 - 13    : 0

Plots

sub_db <- db_hydric %>%
  dplyr::select(chunk_mean_Db_deg_C_diff, CEWL_g_m2h, osmolality_mmol_kg_mean, mass_g, SMI)
PerformanceAnalytics::chart.Correlation(sub_db)
## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

ggplot(db_hydric) +
  geom_point(aes(x = CEWL_g_m2h,
                 y = chunk_mean_Db_deg_C_diff,
                 color = date_chunks)) +
  geom_smooth(aes(x = CEWL_g_m2h,
                 y = chunk_mean_Db_deg_C_diff),
              se = F,
              method = "lm",
              formula = y~x)

ggplot(db_hydric) +
  geom_point(aes(x = osmolality_mmol_kg_mean,
                 y = chunk_mean_Db_deg_C_diff,
                 color = date_chunks)) +
  geom_smooth(aes(x = osmolality_mmol_kg_mean,
                 y = chunk_mean_Db_deg_C_diff),
              se = F,
              method = "lm",
              formula = y~x)
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

Models

# CEWL
db_CEWL_mod <- lmerTest::lmer(data = db_hydric,
                                chunk_mean_Db_deg_C_diff ~ CEWL_g_m2h +
                                  (1|individual_ID))
summary(db_CEWL_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ CEWL_g_m2h + (1 | individual_ID)
##    Data: db_hydric
## 
## REML criterion at convergence: 193.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.62001 -0.24244 -0.01936  0.31381  1.60181 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 5.71708  2.3910  
##  Residual                  0.08465  0.2909  
## Number of obs: 55, groups:  individual_ID, 37
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  2.10721    0.42177 44.35973   4.996  9.6e-06 ***
## CEWL_g_m2h  -0.02517    0.01539 17.28891  -1.635     0.12    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## CEWL_g_m2h -0.349
plot(db_CEWL_mod)

# osmolality
db_osml_mod <- lmerTest::lmer(data = db_hydric,
                                chunk_mean_Db_deg_C_diff ~ osmolality_mmol_kg_mean +
                                  (1|individual_ID))
summary(db_osml_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## chunk_mean_Db_deg_C_diff ~ osmolality_mmol_kg_mean + (1 | individual_ID)
##    Data: db_hydric
## 
## REML criterion at convergence: 189.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.20314 -0.29959 -0.02697  0.34570  1.15461 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 5.689    2.3852  
##  Residual                  0.072    0.2683  
## Number of obs: 54, groups:  individual_ID, 36
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)             -0.634993   1.016775 24.263771  -0.625   0.5381  
## osmolality_mmol_kg_mean  0.006915   0.002553 17.693153   2.708   0.0146 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.920
plot(db_osml_mod)

# mass
db_mass_mod <- lmerTest::lmer(data = db_hydric,
                                chunk_mean_Db_deg_C_diff ~ mass_g +
                                  (1|individual_ID))
summary(db_mass_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ mass_g + (1 | individual_ID)
##    Data: db_hydric
## 
## REML criterion at convergence: 192.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57300 -0.13496 -0.02122  0.17189  1.58553 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 5.73934  2.3957  
##  Residual                  0.08222  0.2867  
## Number of obs: 55, groups:  individual_ID, 37
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)  0.42166    0.90373 33.57494   0.467   0.6438  
## mass_g       0.03762    0.02115 23.77678   1.779   0.0881 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## mass_g -0.899
plot(db_mass_mod)

# SMI
db_SMI_mod <- lmerTest::lmer(data = db_hydric,
                                chunk_mean_Db_deg_C_diff ~ SMI +
                                  (1|individual_ID))
summary(db_SMI_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ SMI + (1 | individual_ID)
##    Data: db_hydric
## 
## REML criterion at convergence: 194.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5647 -0.1435 -0.0334  0.1880  1.5731 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 5.63424  2.3737  
##  Residual                  0.09294  0.3049  
## Number of obs: 55, groups:  individual_ID, 37
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)  1.21859    0.69212 33.91620   1.761   0.0873 .
## SMI          0.01846    0.01624 17.27587   1.137   0.2712  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## SMI -0.824
plot(db_SMI_mod)

# supplemental water tmt
db_tmt_mod <- lmerTest::lmer(data = db_hydric,
                                chunk_mean_Db_deg_C_diff ~ tmt +
                                  (1|individual_ID))
summary(db_tmt_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: chunk_mean_Db_deg_C_diff ~ tmt + (1 | individual_ID)
##    Data: db_hydric
## 
## REML criterion at convergence: 187.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.45783 -0.18050 -0.03887  0.31549  1.48654 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 5.55610  2.3571  
##  Residual                  0.09779  0.3127  
## Number of obs: 55, groups:  individual_ID, 37
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.1925     0.5310 34.8339   4.129 0.000216 ***
## tmtSham Tmt  -0.7086     0.7826 34.6996  -0.905 0.371490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tmtSham Tmt -0.678

Microhabitat ~ Hydric

Prep Data

I actually need to prep the microhabitat data- the metric I will investigate is proportion of time aboveground.

prop_above_hydric <- prop_above_by_indiv %>%
  dplyr::filter(date_chunks %in% c("Apr 26 - May 6", "May 10 - 20"))
micro_hydric <- liz_dat_apr_may_only %>%
  left_join(prop_above_hydric, by = c('date_chunks', 
                              'radio_collar_serial' = 'Serial')) %>%
  # remove NAs
  dplyr::filter(complete.cases(prop_above)) 

summary(micro_hydric)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 03:43:00.00   Min.   :2021-04-23   229-060: 2         
##  1st Qu.:2021-04-23 06:44:00.00   1st Qu.:2021-04-23   252-884: 2         
##  Median :2021-04-24 03:49:00.00   Median :2021-04-24   252-887: 2         
##  Mean   :2021-04-28 00:13:02.60   Mean   :2021-04-28   229-059: 1         
##  3rd Qu.:2021-05-07 04:15:30.00   3rd Qu.:2021-05-08   229-066: 1         
##  Max.   :2021-05-08 07:19:00.00   Max.   :2021-05-08   229-070: 1         
##  NA's   :2                                             (Other):16         
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:25          F-11   : 2    Female:15   Min.   :26.00   Min.   : 85.0  
##  Class :character   F-17   : 2    Male  :10   1st Qu.:32.00   1st Qu.: 94.0  
##  Mode  :character   M-04   : 2                Median :35.60   Median : 97.0  
##                     F-01   : 1                Mean   :36.81   Mean   : 98.4  
##                     F-03   : 1                3rd Qu.:38.00   3rd Qu.:104.0  
##                     F-05   : 1                Max.   :52.70   Max.   :114.0  
##                     (Other):16                                               
##         tmt     tmt_mass_change_g capture_to_msmt  hematocrit_percent
##  Water Tmt:11   Min.   :-3.000    Min.   :  25.0   Min.   :23.00     
##  Sham Tmt :14   1st Qu.:-1.000    1st Qu.:  67.5   1st Qu.:27.75     
##  No Tmt   : 0   Median : 0.000    Median : 154.0   Median :32.50     
##                 Mean   :-0.464    Mean   : 437.7   Mean   :34.46     
##                 3rd Qu.: 0.000    3rd Qu.: 267.5   3rd Qu.:40.25     
##                 Max.   : 1.800    Max.   :1662.0   Max.   :58.00     
##                                   NA's   :2        NA's   :1         
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :24.00   Min.   : 0.650   Min.   :19.07   Min.   :13.12  
##  1st Qu.:32.00   1st Qu.: 8.457   1st Qu.:28.46   1st Qu.:14.63  
##  Median :33.00   Median : 9.495   Median :30.68   Median :16.30  
##  Mean   :32.68   Mean   :10.365   Mean   :29.10   Mean   :19.49  
##  3rd Qu.:34.00   3rd Qu.:12.213   3rd Qu.:31.80   3rd Qu.:22.15  
##  Max.   :36.00   Max.   :21.673   Max.   :33.07   Max.   :35.83  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.205   Min.   :1.415   Min.   :316.0           April:16   1 :16  
##  1st Qu.:3.881   1st Qu.:2.990   1st Qu.:340.2           May  : 9   3 : 9  
##  Median :4.411   Median :3.703   Median :357.8           July : 0   12: 0  
##  Mean   :4.109   Mean   :3.358   Mean   :359.4                             
##  3rd Qu.:4.701   3rd Qu.:4.007   3rd Qu.:378.6                             
##  Max.   :5.049   Max.   :4.309   Max.   :402.0                             
##                                  NA's   :1                                 
##     week_num    ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1.00   Min.   :2021-04-24   Gravid    : 9   Min.   :3.000  
##  1st Qu.:1.00   1st Qu.:2021-04-24   Not Gravid: 5   1st Qu.:3.000  
##  Median :1.00   Median :2021-05-01   NA's      :11   Median :4.000  
##  Mean   :1.72   Mean   :2021-05-01                   Mean   :3.667  
##  3rd Qu.:3.00   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3.00   Max.   :2021-05-08                   Max.   :4.000  
##                 NA's   :11                           NA's   :16     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage   
##  Min.   :0.640           Min.   :1.82                 small round: 1  
##  1st Qu.:0.870           1st Qu.:2.45                 large round: 6  
##  Median :0.990           Median :2.88                 soft       : 1  
##  Mean   :1.043           Mean   :2.95                 soft oblong: 0  
##  3rd Qu.:1.200           3rd Qu.:3.48                 firm oblong: 1  
##  Max.   :1.510           Max.   :4.07                 hard oblong: 0  
##  NA's   :16              NA's   :16                   NA's       :16  
##    dev_point          SMI                date_chunks          microhabitat
##  Min.   :1.000   Min.   :28.35   Apr 26 - May 6:16   Open           : 0   
##  1st Qu.:2.000   1st Qu.:32.48   May 10 - 20   : 9   Shade - Partial: 0   
##  Median :2.000   Median :34.49   May 21 - 31   : 0   Shade - Full   : 0   
##  Mean   :2.333   Mean   :35.16   Jun 1 - 11    : 0   Burrow         :25   
##  3rd Qu.:2.000   3rd Qu.:36.44   Jun 12 - 22   : 0                        
##  Max.   :5.000   Max.   :45.63   Jun 23 - Jul 3: 0                        
##  NA's   :16                      Jul 4 - 13    : 0                        
##   n_obs_micro    n_obs_total  MH_use_proportion   prop_above    
##  Min.   :1.00   Min.   :4.0   Min.   :0.1429    Min.   :0.2500  
##  1st Qu.:1.00   1st Qu.:4.0   1st Qu.:0.2000    1st Qu.:0.6000  
##  Median :1.00   Median :5.0   Median :0.2500    Median :0.7500  
##  Mean   :1.68   Mean   :5.2   Mean   :0.3287    Mean   :0.6713  
##  3rd Qu.:2.00   3rd Qu.:6.0   3rd Qu.:0.4000    3rd Qu.:0.8000  
##  Max.   :4.00   Max.   :7.0   Max.   :0.7500    Max.   :0.8571  
## 

Plots

sub_micro <- micro_hydric %>%
  dplyr::select(prop_above, CEWL_g_m2h, osmolality_mmol_kg_mean, mass_g, SMI)
PerformanceAnalytics::chart.Correlation(sub_micro)
## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

ggplot(MH_use_by_indiv) + 
    geom_bar(aes(x = date_chunks,
                 y = MH_use_proportion,
                 fill = microhabitat
                 ),
             position = "fill", 
             stat = "identity") +
  facet_wrap(~Serial) +
  theme_classic()

ggplot(micro_hydric) +
  geom_point(aes(x = CEWL_g_m2h,
                 y = prop_above,
                 color = date_chunks)) +
  geom_smooth(aes(x = CEWL_g_m2h,
                 y = prop_above),
              se = F,
              method = "lm",
              formula = y~x)

ggplot(micro_hydric) +
  geom_point(aes(x = osmolality_mmol_kg_mean,
                 y = prop_above,
                 color = date_chunks)) +
  geom_smooth(aes(x = osmolality_mmol_kg_mean,
                 y = prop_above),
              se = F,
              method = "lm",
              formula = y~x)
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

Models

# CEWL
micro_CEWL_mod <- lmerTest::lmer(data = micro_hydric,
                                prop_above ~ CEWL_g_m2h +
                                  (1|individual_ID))
summary(micro_CEWL_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ CEWL_g_m2h + (1 | individual_ID)
##    Data: micro_hydric
## 
## REML criterion at convergence: -9.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.99464 -0.18300  0.07242  0.16112  0.99194 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.03197  0.17881 
##  Residual                  0.00143  0.03781 
## Number of obs: 25, groups:  individual_ID, 22
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  0.791099   0.059783  7.601336   13.23 1.61e-06 ***
## CEWL_g_m2h  -0.012504   0.004465  2.748072   -2.80    0.075 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## CEWL_g_m2h -0.759
plot(micro_CEWL_mod)

# osmolality
micro_osml_mod <- lmerTest::lmer(data = micro_hydric,
                                prop_above ~ osmolality_mmol_kg_mean +
                                  (1|individual_ID))
summary(micro_osml_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ osmolality_mmol_kg_mean + (1 | individual_ID)
##    Data: micro_hydric
## 
## REML criterion at convergence: -3.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0054 -0.4661  0.1346  0.3529  0.9514 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.026216 0.16191 
##  Residual                  0.005785 0.07606 
## Number of obs: 24, groups:  individual_ID, 21
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)              1.031922   0.420963  8.191928   2.451   0.0392 *
## osmolality_mmol_kg_mean -0.001043   0.001166  8.075337  -0.895   0.3966  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.996
plot(micro_osml_mod)

# mass
micro_mass_mod <- lmerTest::lmer(data = micro_hydric,
                                prop_above ~ mass_g +
                                  (1|individual_ID))
summary(micro_mass_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ mass_g + (1 | individual_ID)
##    Data: micro_hydric
## 
## REML criterion at convergence: -6.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.98776 -0.41604  0.06801  0.40571  0.95433 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.025796 0.16061 
##  Residual                  0.006184 0.07864 
## Number of obs: 25, groups:  individual_ID, 22
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)  0.524904   0.187152 21.190865   2.805   0.0106 *
## mass_g       0.003797   0.004978 21.259937   0.763   0.4540  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## mass_g -0.979
plot(micro_mass_mod)

# SMI
micro_SMI_mod <- lmerTest::lmer(data = micro_hydric,
                                prop_above ~ SMI +
                                  (1|individual_ID))
summary(micro_SMI_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ SMI + (1 | individual_ID)
##    Data: micro_hydric
## 
## REML criterion at convergence: -7.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0319 -0.3658  0.1095  0.3462  0.9900 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.028661 0.16930 
##  Residual                  0.004752 0.06893 
## Number of obs: 25, groups:  individual_ID, 22
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)  0.472810   0.306820 22.773223   1.541    0.137
## SMI          0.005427   0.008619 22.720145   0.630    0.535
## 
## Correlation of Fixed Effects:
##     (Intr)
## SMI -0.992
plot(micro_SMI_mod)

# supplemental hydration tmt
micro_tmt_mod <- lmerTest::lmer(data = micro_hydric,
                                prop_above ~ tmt +
                                  (1|individual_ID))
summary(micro_tmt_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop_above ~ tmt + (1 | individual_ID)
##    Data: micro_hydric
## 
## REML criterion at convergence: -12.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0381 -0.2972  0.1325  0.3188  0.9215 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 0.027290 0.1652  
##  Residual                  0.005315 0.0729  
## Number of obs: 25, groups:  individual_ID, 22
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.63293    0.05685 19.84699  11.133 5.52e-10 ***
## tmtSham Tmt  0.05782    0.07687 19.72667   0.752    0.461    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tmtSham Tmt -0.740

All Anovas

Just so they’re all in one place, get all the anovas from the models above to run:

CEWL

anova(max_temp_CEWL_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##             Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## CEWL_g_m2h 0.24422 0.24422     1 42.422  0.3076 0.5821
anova(db_CEWL_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##             Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## CEWL_g_m2h 0.22482 0.22482     1 17.734  2.6559 0.1208
anova(micro_CEWL_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##              Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## CEWL_g_m2h 0.008691 0.008691     1 2.8375  6.0795 0.09528 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Osmolality

anova(max_temp_osml_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## osmolality_mmol_kg_mean  2.073   2.073     1 46.788  2.9908 0.09033 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(max_temp_osml_mod)
anova(db_osml_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                          Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## osmolality_mmol_kg_mean 0.52367 0.52367     1 17.855  7.2737 0.01482 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(db_osml_mod)
anova(micro_osml_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                            Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 0.0030997 0.0030997     1 7.7022  0.5358 0.4858

Body Mass

anova(max_temp_mass_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##        Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## mass_g 1.1961  1.1961     1 41.917  1.5485 0.2203
anova(db_mass_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##         Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## mass_g 0.24808 0.24808     1 24.158  3.0171 0.09513 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(micro_mass_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##           Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## mass_g 0.0034519 0.0034519     1 21.557  0.5581 0.4631

Body Condition

anova(max_temp_SMI_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##     Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## SMI 3.3291  3.3291     1 41.314  4.3023 0.04433 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(max_temp_SMI_mod)
anova(db_SMI_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##      Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## SMI 0.11927 0.11927     1 17.676  1.2833 0.2724
anova(micro_SMI_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##        Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## SMI 0.0015565 0.0015565     1 22.753  0.3276 0.5727

Correlation “Matrix”

For each model, get the marginal Rsq for the amount of variance that’s explained by fixed effects only. Then, we will use that to calculate R, and make a figure. I put “matrix” in quotations because it will be a subset of what would be presented for a full correlation matrix.

corrs <- data.frame(hydric_var = c(rep("CEWL", 3), rep("Plasma Osmolality", 3), 
                                   #rep("Body Mass", 3), 
                                   rep("Body Condition", 3)),
                    thermal_var = c(rep(c("Max Temp", "Db", "Micro Use"), 3)),
                    R = c(
                      # CEWL
                      cor(max_temps_hydric$CEWL_g_m2h,
                          max_temps_hydric$chunk_avg_Tb_percentile_80_90),
                      cor(db_hydric$CEWL_g_m2h, db_hydric$chunk_mean_Db_deg_C_diff),
                      cor(micro_hydric$CEWL_g_m2h, micro_hydric$prop_above),
                      # osmolality
                      cor(max_temps_hydric$osmolality_mmol_kg_mean,
                          max_temps_hydric$chunk_avg_Tb_percentile_80_90,
                          use = "complete.obs"),
                      cor(db_hydric$osmolality_mmol_kg_mean,
                          db_hydric$chunk_mean_Db_deg_C_diff,
                          use = "complete.obs"),
                      cor(micro_hydric$osmolality_mmol_kg_mean, micro_hydric$prop_above,
                          use = "complete.obs"),
                      # mass
                      # cor(max_temps_hydric$mass_g,
                      #     max_temps_hydric$chunk_avg_Tb_percentile_80_90),
                      # cor(db_hydric$mass_g, db_hydric$chunk_mean_Db_deg_C_diff),
                      # cor(micro_hydric$mass_g, micro_hydric$prop_above),
                      # body condition
                      cor(max_temps_hydric$SMI,
                          max_temps_hydric$chunk_avg_Tb_percentile_80_90),
                      cor(db_hydric$SMI, db_hydric$chunk_mean_Db_deg_C_diff),
                      cor(micro_hydric$SMI, micro_hydric$prop_above))) %>%
  mutate(hydric_var = factor(hydric_var, levels = c("CEWL", "Plasma Osmolality",
                                                    #"Body Mass", 
                                                    "Body Condition")))
                   
corrs
##          hydric_var thermal_var           R
## 1              CEWL    Max Temp  0.03561470
## 2              CEWL          Db  0.24284741
## 3              CEWL   Micro Use -0.11107187
## 4 Plasma Osmolality    Max Temp  0.16506018
## 5 Plasma Osmolality          Db  0.03476896
## 6 Plasma Osmolality   Micro Use -0.14174829
## 7    Body Condition    Max Temp  0.26269876
## 8    Body Condition          Db -0.21638149
## 9    Body Condition   Micro Use  0.03188603

Plot it:

removed body mass because body condition was more relevant

ggplot(corrs) +
  aes(x = hydric_var,
      y = thermal_var,
      fill = R) +
  geom_tile() +
  scale_fill_distiller(palette = "RdYlBu", limits = c(-1, 1),
                       name = "Correlation\nCoefficient") +
  savs_ggplot_theme2 +
  scale_y_discrete(name = NULL, labels = c("Thermoregulatory\n            Accuracy",
                                           bquote("Maximum T"[b]),
                                           "Microhabitat Use"
                                           )) +
  scale_x_discrete(name = NULL, labels = c("CEWL", "Plasma\nOsmolality",
                                           #"Body\nMass", 
                                           "Body\nCondition")) +
  # CEWL
  annotate("text", x = 1, y = 1, label = "-0.03", 
           size = 2, family = "sans", color = "black") +
  annotate("text", x = 1, y = 2, label = "-0.02", 
           size = 2, family = "sans", color = "black") +
  annotate("text", x = 1, y = 3, label = "-0.01", 
           size = 2, family = "sans", color = "black") +
  # plasma osmolality
  annotate("text", x = 2, y = 1, label = "0.007*", 
           size = 3, family = "sans", color = "black") +
  annotate("text", x = 2, y = 2, label = "0.01", 
           size = 2, family = "sans", color = "black") +
  annotate("text", x = 2, y = 3, label = "-0.001", 
           size = 2, family = "sans", color = "black") +
  # body mass
  # annotate("text", x = 3, y = 1, label = "0.04", 
  #          size = 2, family = "sans", color = "black") +
  # annotate("text", x = 3, y = 2, label = "0.03", 
  #          size = 2, family = "sans", color = "black") +
  # annotate("text", x = 3, y = 3, label = "0.004", 
  #          size = 2, family = "sans", color = "black") +
  # body condition
  annotate("text", x = 3, y = 1, label = "0.02", 
           size = 2, family = "sans", color = "black") +
  annotate("text", x = 3, y = 2, label = "0.08*", 
           size = 3, family = "sans", color = "black") +
  annotate("text", x = 3, y = 3, label = "0.005", 
           size = 2, family = "sans", color = "black") -> corr_plot

# arrange to get legend centered
ggarrange(corr_plot,
          ncol = 1, nrow = 1,
          common.legend = TRUE,
          legend = "bottom"
          ) -> corr_plot_arrange

# export
ggsave(filename = "thermal_hydric_corr_plot.pdf",
       plot = corr_plot_arrange,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 100, height = 80)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## metrics unknown for character 0xa
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font metrics unknown for character 0xa

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font metrics unknown for character 0xa

Hydric Conclusion

There are basically no detectable relationships between water balance (CEWL, osml, or mass) and thermoregulation or microhabitat use in BNLL.

Microhabitat Use ONLY

Prep Data

For logistic regression, I need to have the full, raw dataset (vs proportions by group for plots), so I will actually use the ‘microhab’ dataframe, created in the “Load Data” > “Microhabitat Use” section.

Basic Plots

ggplot(MH_use_across_indivs_sex) + 
    geom_bar(aes(x = date_chunks,
                 y = MH_use_proportion,
                 fill = microhabitat
                 ),
             position = "fill", 
             stat = "identity") +
  theme_classic() +
  facet_wrap(~sex_M_F)

ggplot(MH_use_across_indivs) + 
    geom_bar(aes(x = date_chunks,
                 y = MH_use_proportion,
                 fill = microhabitat
                 ),
             position = "fill", 
             stat = "identity") +
  theme_classic()

Logistic Regression

Start with only ~ time period/chunk:

micro_glmm <- lme4::glmer(data = microhab, 
                           above_below ~ date_chunks +
                            (1|individual_ID),
                           family="binomial")
summary(micro_glmm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: above_below ~ date_chunks + (1 | individual_ID)
##    Data: microhab
## 
##      AIC      BIC   logLik deviance df.resid 
##   1688.7   1732.6   -836.4   1672.7     1767 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1617 -0.6023 -0.1103  0.4143 11.2158 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 2.023    1.422   
## Number of obs: 1775, groups:  individual_ID, 39
## 
## Fixed effects:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.709551   0.315963  -5.411 6.28e-08 ***
## date_chunksMay 10 - 20     0.419879   0.301382   1.393    0.164    
## date_chunksMay 21 - 31    -0.004787   0.315245  -0.015    0.988    
## date_chunksJun 1 - 11      0.429209   0.316201   1.357    0.175    
## date_chunksJun 12 - 22     2.669814   0.307605   8.679  < 2e-16 ***
## date_chunksJun 23 - Jul 3  3.159663   0.291237  10.849  < 2e-16 ***
## date_chunksJul 4 - 13      4.177164   0.335053  12.467  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) d_M1-2 d_M2-3 d_J1-1 d_J1-2 d_2-J3
## dt_chM10-20 -0.484                                   
## dt_chM21-31 -0.470  0.550                            
## dt_chnJ1-11 -0.477  0.537  0.560                     
## dt_chJ12-22 -0.516  0.550  0.565  0.642              
## dt_chJ23-J3 -0.550  0.580  0.596  0.678  0.799       
## dt_chnJ4-13 -0.490  0.505  0.519  0.592  0.715  0.784
pairs(emmeans(micro_glmm, "date_chunks"))
##  contrast                            estimate    SE  df z.ratio p.value
##  (Apr 26 - May 6) - (May 10 - 20)    -0.41988 0.301 Inf  -1.393  0.8058
##  (Apr 26 - May 6) - (May 21 - 31)     0.00479 0.315 Inf   0.015  1.0000
##  (Apr 26 - May 6) - (Jun 1 - 11)     -0.42921 0.316 Inf  -1.357  0.8244
##  (Apr 26 - May 6) - (Jun 12 - 22)    -2.66981 0.308 Inf  -8.679  <.0001
##  (Apr 26 - May 6) - (Jun 23 - Jul 3) -3.15966 0.291 Inf -10.849  <.0001
##  (Apr 26 - May 6) - (Jul 4 - 13)     -4.17716 0.335 Inf -12.467  <.0001
##  (May 10 - 20) - (May 21 - 31)        0.42467 0.293 Inf   1.450  0.7743
##  (May 10 - 20) - (Jun 1 - 11)        -0.00933 0.298 Inf  -0.031  1.0000
##  (May 10 - 20) - (Jun 12 - 22)       -2.24993 0.289 Inf  -7.784  <.0001
##  (May 10 - 20) - (Jun 23 - Jul 3)    -2.73978 0.272 Inf -10.081  <.0001
##  (May 10 - 20) - (Jul 4 - 13)        -3.75728 0.318 Inf -11.813  <.0001
##  (May 21 - 31) - (Jun 1 - 11)        -0.43400 0.296 Inf  -1.465  0.7658
##  (May 21 - 31) - (Jun 12 - 22)       -2.67460 0.291 Inf  -9.206  <.0001
##  (May 21 - 31) - (Jun 23 - Jul 3)    -3.16445 0.273 Inf -11.577  <.0001
##  (May 21 - 31) - (Jul 4 - 13)        -4.18195 0.319 Inf -13.090  <.0001
##  (Jun 1 - 11) - (Jun 12 - 22)        -2.24060 0.264 Inf  -8.488  <.0001
##  (Jun 1 - 11) - (Jun 23 - Jul 3)     -2.73045 0.245 Inf -11.156  <.0001
##  (Jun 1 - 11) - (Jul 4 - 13)         -3.74796 0.295 Inf -12.723  <.0001
##  (Jun 12 - 22) - (Jun 23 - Jul 3)    -0.48985 0.191 Inf  -2.569  0.1356
##  (Jun 12 - 22) - (Jul 4 - 13)        -1.50735 0.244 Inf  -6.184  <.0001
##  (Jun 23 - Jul 3) - (Jul 4 - 13)     -1.01750 0.210 Inf  -4.848  <.0001
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 7 estimates

pairwise comparison groups for when I make figures:

  • Apr 26 - May 6 = A
  • May 10 - 20 = A
  • May 21 - 31 = A
  • Jun 1 - 11 = A
  • Jun 12 - 22 = B
  • Jun 23 - Jul 3 = B
  • Jul 4 - 13 = C (wow!!)

Compare a model with only ~ chunks vs chunks*sex:

micro_glmm_sex <- lme4::glmer(data = microhab, 
                           above_below ~ date_chunks*sex_M_F +
                            (1|individual_ID),
                           family="binomial")
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0650045 (tol = 0.002, component 1)
summary(micro_glmm_sex)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: above_below ~ date_chunks * sex_M_F + (1 | individual_ID)
##    Data: microhab
## 
##      AIC      BIC   logLik deviance df.resid 
##   1677.1   1759.4   -823.6   1647.1     1760 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0934 -0.5633 -0.0851  0.4084 14.0614 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 1.828    1.352   
## Number of obs: 1775, groups:  individual_ID, 39
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            -1.4451     0.4263  -3.390 0.000699 ***
## date_chunksMay 10 - 20                  1.1237     0.4159   2.702 0.006899 ** 
## date_chunksMay 21 - 31                  0.6417     0.4552   1.410 0.158613    
## date_chunksJun 1 - 11                   1.3148     0.4952   2.655 0.007936 ** 
## date_chunksJun 12 - 22                  2.8696     0.4985   5.756 8.61e-09 ***
## date_chunksJun 23 - Jul 3               3.4974     0.4714   7.420 1.17e-13 ***
## date_chunksJul 4 - 13                   3.4882     0.5530   6.307 2.84e-10 ***
## sex_M_FMale                            -0.4569     0.6125  -0.746 0.455696    
## date_chunksMay 10 - 20:sex_M_FMale     -1.5098     0.6171  -2.447 0.014416 *  
## date_chunksMay 21 - 31:sex_M_FMale     -1.2489     0.6350  -1.967 0.049199 *  
## date_chunksJun 1 - 11:sex_M_FMale      -1.5622     0.6467  -2.416 0.015706 *  
## date_chunksJun 12 - 22:sex_M_FMale     -0.4402     0.6282  -0.701 0.483534    
## date_chunksJun 23 - Jul 3:sex_M_FMale  -0.6205     0.5934  -1.046 0.295746    
## date_chunksJul 4 - 13:sex_M_FMale       0.6935     0.6864   1.010 0.312392    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0650045 (tol = 0.002, component 1)
car::Anova(micro_glmm_sex, type = 2, 
           test.statistic = "Chisq", error.estimate = "pearson")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: above_below
##                        Chisq Df Pr(>Chisq)    
## date_chunks         283.6728  6  < 2.2e-16 ***
## sex_M_F               5.3173  1   0.021115 *  
## date_chunks:sex_M_F  20.3379  6   0.002411 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(micro_glmm, micro_glmm_sex)
## Data: microhab
## Models:
## micro_glmm: above_below ~ date_chunks + (1 | individual_ID)
## micro_glmm_sex: above_below ~ date_chunks * sex_M_F + (1 | individual_ID)
##                npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## micro_glmm        8 1688.7 1732.6 -836.36   1672.7                         
## micro_glmm_sex   15 1677.2 1759.4 -823.57   1647.2 25.581  7  0.0005982 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Okay, adding sex was technically a better model based on some metrics, but how many time periods was microhabitat use actually different between the sexes?

emmeans_micro_sex_glmm_sex_diffs <- data.frame(emmeans(micro_glmm_sex, 
                          pairwise ~ sex_M_F | date_chunks)$contrasts) %>%
  dplyr::select(date_chunks, contrast, p.value)
emmeans_micro_sex_glmm_sex_diffs
##      date_chunks      contrast     p.value
## 1 Apr 26 - May 6 Female - Male 0.455695676
## 2    May 10 - 20 Female - Male 0.001642991
## 3    May 21 - 31 Female - Male 0.007776951
## 4     Jun 1 - 11 Female - Male 0.001975646
## 5    Jun 12 - 22 Female - Male 0.150030024
## 6 Jun 23 - Jul 3 Female - Male 0.066251295
## 7     Jul 4 - 13 Female - Male 0.726818004

May - June, microhabitat use was significantly different based on sex. Females were spending more time in burrows than males (based on the plots above).

To keep it simple, let’s keep the main figure to showing both males and females aggregated over time, but I can also make a secondary figure showing their differences in microhabitat use specifically in May-early June.

Probabilities from GLMM

Now, I want to use the model estimates/coefficients to calculate the predicted probability of lizards being belowground. First, save the estimates and standard errors. I can do this by getting the model-estimated means. Then, I can use this equation to calculate the model-predicted probability of being belowground:

p = exp(emmean) / (1 + exp(emmean))

# one set for M-F differences
emmeans_probabilities_MF <- data.frame(emmeans(micro_glmm_sex, 
                          pairwise ~ sex_M_F | date_chunks)$emmeans) %>% 
  # calculate prob
  mutate(exp_emmean = exp(emmean),
         prob = exp_emmean / (1 + exp_emmean),
         # and confidence intervals, in case helpful
         exp_LCL = exp(asymp.LCL),
         prob_LCL = exp_LCL / (1 + exp_LCL),
         exp_UCL = exp(asymp.UCL),
         prob_UCL = exp_UCL / (1 + exp_UCL)
         ) %>% 
  # only save the probabilities for when M & F were sig diff
  dplyr::filter(date_chunks %in% c("May 10 - 20", "May 21 - 31", "Jun 1 - 11"))
emmeans_probabilities_MF
##   sex_M_F date_chunks     emmean        SE  df asymp.LCL   asymp.UCL exp_emmean
## 1  Female May 10 - 20 -0.3214039 0.4264786 Inf -1.157286  0.51447873 0.72513033
## 2    Male May 10 - 20 -2.2880805 0.4562992 Inf -3.182410 -1.39375044 0.10146103
## 3  Female May 21 - 31 -0.8034043 0.4608742 Inf -1.706701  0.09989252 0.44780191
## 4    Male May 21 - 31 -2.5091531 0.4450839 Inf -3.381502 -1.63680471 0.08133709
## 5  Female  Jun 1 - 11 -0.1303629 0.4964374 Inf -1.103362  0.84263653 0.87777684
## 6    Male  Jun 1 - 11 -2.1494058 0.4224566 Inf -2.977406 -1.32140601 0.11655339
##         prob    exp_LCL   prob_LCL   exp_UCL  prob_UCL
## 1 0.42033365 0.31433798 0.23916069 1.6727663 0.6258558
## 2 0.09211495 0.04148553 0.03983304 0.2481429 0.1988097
## 3 0.30929778 0.18146343 0.15359208 1.1050521 0.5249524
## 4 0.07521900 0.03399637 0.03287862 0.1946009 0.1629003
## 5 0.46745536 0.33175375 0.24911043 2.3224822 0.6990202
## 6 0.10438676 0.05092478 0.04845712 0.2667600 0.2105845
# one set for overall probabilities
emmeans_probabilities <- data.frame(emmeans(micro_glmm, "date_chunks")) %>% 
  # calculate prob
  mutate(exp_emmean = exp(emmean),
         prob = exp_emmean / (1 + exp_emmean),
         # and confidence intervals, in case helpful
         exp_LCL = exp(asymp.LCL),
         prob_LCL = exp_LCL / (1 + exp_LCL),
         exp_UCL = exp(asymp.UCL),
         prob_UCL = exp_UCL / (1 + exp_UCL)
         ) %>% 
  # only save the probabilities for when M & F were the same
  dplyr::filter(date_chunks %nin% c("May 10 - 20", "May 21 - 31", "Jun 1 - 11")
                ) %>% 
  # add "sex" col
  mutate(sex_M_F = "Both")
emmeans_probabilities
##      date_chunks     emmean        SE  df  asymp.LCL asymp.UCL exp_emmean
## 1 Apr 26 - May 6 -1.7095505 0.3159632 Inf -2.3288271 -1.090274  0.1809471
## 2    Jun 12 - 22  0.9602636 0.3067858 Inf  0.3589745  1.561553  2.6123850
## 3 Jun 23 - Jul 3  1.4501123 0.2887028 Inf  0.8842652  2.015959  4.2635933
## 4     Jul 4 - 13  2.4676138 0.3292154 Inf  1.8223634  3.112864 11.7942692
##        prob    exp_LCL   prob_LCL    exp_UCL  prob_UCL sex_M_F
## 1 0.1532220 0.09740994 0.08876349  0.3361244 0.2515667    Both
## 2 0.7231746 1.43186024 0.58879216  4.7662161 0.8265760    Both
## 3 0.8100157 2.42120453 0.70770529  7.5079274 0.8824626    Both
## 4 0.9218400 6.18646237 0.86084948 22.4853524 0.9574203    Both
# put them together
barplot_probs <- emmeans_probabilities %>% 
  rbind(emmeans_probabilities_MF) %>% 
  dplyr::select(date_chunks, sex_M_F,
                prob, prob_LCL, prob_UCL) %>% 
  mutate(sex_M_F = factor(sex_M_F))

summary(barplot_probs)
##          date_chunks   sex_M_F       prob            prob_LCL      
##  Apr 26 - May 6:1    Both  :4   Min.   :0.07522   Min.   :0.03288  
##  May 10 - 20   :2    Female:3   1st Qu.:0.11660   1st Qu.:0.05853  
##  May 21 - 31   :2    Male  :3   Median :0.36482   Median :0.19638  
##  Jun 1 - 11    :2               Mean   :0.40771   Mean   :0.30091  
##  Jun 12 - 22   :1               3rd Qu.:0.65924   3rd Qu.:0.50387  
##  Jun 23 - Jul 3:1               Max.   :0.92184   Max.   :0.86085  
##  Jul 4 - 13    :1                                                  
##     prob_UCL     
##  Min.   :0.1629  
##  1st Qu.:0.2208  
##  Median :0.5754  
##  Mean   :0.5340  
##  3rd Qu.:0.7947  
##  Max.   :0.9574  
## 

Pretty Barplot

ggplot() + 
  
  # raw observed proportions
  geom_bar(data = MH_use_across_indivs,
    aes(x = date_chunks,
        y = MH_use_proportion,
        fill = microhabitat
                 ),
             position = "fill", 
             stat = "identity"
    ) +
  
  # model-predicted probabilities
  geom_errorbar(data = barplot_probs,
             aes(
               x = date_chunks,
               ymin = prob_LCL,
               ymax = prob_UCL,
             ),
             color = "black",
             width = 0.1
             ) +
  
  geom_point(data = barplot_probs,
             aes(
               x = date_chunks,
               y = prob,
               shape = sex_M_F,
             ),
             color = "black",
             size = 2
             ) +
  
  # labels
  xlab("Time Period") +
  ylab("Microhabitat Use (proportion of time)") +
  
  # hydric phys msmsts
  geom_vline(xintercept = 0.5, # first weekend of msmts
             size = 0.5, # default size is 0.5
             linetype = "dashed") +
  geom_vline(xintercept = 1.5, # second weekend of msmts
             size = 0.5, 
             linetype = "dashed") +
  geom_vline(xintercept = 7.5, # last few msmts
             size = 0.5, 
             linetype = "dashed") +
  
  # pretties 
  scale_shape_manual(values = c(17, 16, 15), name = "Lizard Sex") +
  scale_fill_manual(values = my_brew_cols_4, name = "Microhabitat",
                    labels = c("Open", "Partial Shade", "Full Shade", "Burrow")) +
  scale_x_discrete(labels = c("Apr 26\n-May 6",
                              "May\n10-20",
                              "May\n21-31",
                              "Jun\n1-11",
                              "Jun\n12-22",
                              "Jun 23\n-Jul 3",
                              "Jul\n4-13"
                              ),
                   name = NULL) +
  scale_y_continuous(limits = c(0, 1.1), expand = c(0,0),
                     breaks = c(seq(0, 1, by = 0.2)),
                     labels = c(seq(0, 1, by = 0.2))) +
  
  # significance annotation for time
  annotate(geom = "text", x = 1, y = 1.08, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 1.08, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 1.08, label = "a", size = 3) +
  annotate(geom = "text", x = 4, y = 1.08, label = "a", size = 3) +
  annotate(geom = "text", x = 5, y = 1.08, label = "b", size = 3) +
  annotate(geom = "text", x = 6, y = 1.08, label = "b", size = 3) +
  annotate(geom = "text", x = 7, y = 1.08, label = "c", size = 3) +
  # significance annotation for sex
  annotate(geom = "text", x = 2, y = 1.015, label = "*", size = 5) +
  annotate(geom = "text", x = 3, y = 1.015, label = "*", size = 5) +
  annotate(geom = "text", x = 4, y = 1.015, label = "*", size = 5) +
  savs_ggplot_theme -> fancy_micro_plot
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
fancy_micro_plot

# arrange to get legend centered
ggarrange(fancy_micro_plot,
          ncol = 1, nrow = 1,
          common.legend = TRUE,
          legend = "right"
          ) -> fancy_micro_plot_arrange

ggsave(filename = "microhabitat_use.pdf",
       plot = fancy_micro_plot_arrange,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 160, height = 80)